Breeding Bird Surveys - Cleveland Metroparks

Nathan W. Byer

5/16/2022

The dataset

Tim Krynak, Erik Shaffer, Jennifer Brumsfield, and other Cleveland Metroparks and affiliated colleagues conducted Breeding Bird Surveys (BBS) at a total of 211 routes between 2017 and 2021. 164 of these routes overlapped with long-term vegetation monitoring efforts through the Cleveland Metroparks Plant Community Assessment Program (PCAP). This provided an excellent opportunity to explore the potential relationships between vegetative composition, landscape structure, and breeding bird ecology and community composition.

Initial processing

First, we just read in the full BBS dataset. For our purposes, we will remove any BBS observations not paired with PCAP plots. We can also visualize the distribution of survey effort through time.

CMP_bbs_pcap_join<-read.csv("CMP_bbs_pcapjoin.csv",header=T)
CMP_bbs_pcap_join_pcaponly<-CMP_bbs_pcap_join[!is.na(CMP_bbs_pcap_join$Plot.ID),]

datetime<-strptime(CMP_bbs_pcap_join_pcaponly$survey_cre,format="%m/%d/%Y %H:%M")
hist(datetime,breaks="months")

Next, we need to start grouping this dataset for downsteram analysis. This involves summarizing counts of each species for each plot. We also replace any NA values in the data frame with zeroes - this will be important for downstream applications, as we want to distinguish between unobserved species at a plot and check (typically coded 0s) and checks that were “skipped” (NAs).

CMP_bbs_pcap_join_pcaponly_pcapsum<-CMP_bbs_pcap_join_pcaponly %>% 
  group_by(simp_plot, reservatio,species_4_) %>%
  summarise(n=n(),longitude=mean(xcoord), latitude=mean(ycoord),vibifq=mean(vibi_fq),vibifqnotrees=mean(vibi_fq_no_trees),carex=mean(carex_metric_value),cyper=mean(cyperaceae_metric_value),
            dicot=mean(dicot_metric_value),shade=mean(shade_metric_value),shrub=mean(shrub_metric_value),
            hydro=mean(hydrophyte_metric_value),svp=mean(svp_metric_value),ap=mean(ap_ratio_metric_value),
            fqi=mean(fqai_score),bryo=mean(bryophyte_metric_value),hydro_rel=mean(hydrophyte_rel_cov_metric_value),
            hydro_rel_notrees=mean(hydrophyte_rel_cov_no_trees),
            sens_rel=mean(sensitive_rel_cov_metric_value),
            sens_rel_notree=mean(sensitive_rel_cov_no_trees),
            tol_rel=mean(tolerant_rel_cov_metric_value),
            tol_rel_notree=mean(tolerant_rel_cov_no_trees),
            inv=mean(invasive_graminoids_metric_value),
            small_tree=mean(small_tree_metric_value),
            subcanopy=mean(subcanopy_iv),
            canopy=mean(canopy_iv.1)
            ) %>%
  spread(species_4_,n) %>%
  replace(is.na(.),0)

CMP_bbs_pcap_join_pcaponly_pcapsum<-column_to_rownames(CMP_bbs_pcap_join_pcaponly_pcapsum,var="simp_plot")

We can also take a look at this dataset to ensure that it is what we expect.

head(CMP_bbs_pcap_join_pcaponly_pcapsum)
##      reservatio longitude latitude   vibifq vibifqnotrees carex cyper dicot shade shrub hydro svp         ap      fqi        bryo  hydro_rel hydro_rel_notrees   sens_rel sens_rel_notree    tol_rel tol_rel_notree inv small_tree  subcanopy    canopy V1 ACFL ALFL AMCR AMGO AMKE AMRE AMRO AMWO BADO
## 1002         HI   2181968 561930.7 47.26427      24.35064     9     9    30    25     1     9   1 0.05882353 19.47063 0.005272844 0.01581853        0.03168627 0.03962464     0.000313725 0.41913886     0.76606536   0 0.13936430 0.03419060 0.1922540  0    1    0    1    2    0    2    0    0    0
## 1003         MS   2150599 613095.3 76.16889      52.67227     4     4    22    22     1     2   1 0.13333333 22.98097 0.009033832 0.01806766        0.06677053 0.48845928     0.051747162 0.11096557     0.04284442   0 0.04845815 0.03497367 0.2054565  1    5    0    0    4    0    0    4    0    0
## 1004         BE   2227458 626184.1 37.88422      35.10323     2     2    21    23     1     9   4 0.19047619 15.08659 0.039605806 0.42634485        0.55871039 0.01397852     0.009159187 0.19872796     0.26042621   0 0.10869565 0.46386133 0.2174377  1    5    0    3    1    0    0    5    0    0
## 1005         HI   2191849 580242.5 82.91987      65.84743    11    11    47    52     2    18   3 0.08000000 30.51229 0.007461016 0.06013579        0.12854181 0.16854436     0.088512041 0.13320401     0.19685290   0 0.11436170 0.10777783 0.1757489  0    3    0    2    2    0    0    2    0    0
## 1006         BW   2116879 635151.5 72.07176      64.74055     3     3    21    23     1     6   2 0.06666667 26.58228 0.008306874 0.11007532        0.25470932 0.25985749     0.046004015 0.12386473     0.06877109   0 0.13235294 0.22866471 0.1707812  0    2    0    1    1    0    0    4    0    0
## 1008         WC   2190113 621267.9 56.30316      25.13574     4     4    11    13     0     3   0 0.25000000 16.84286 0.010954302 0.05477151        0.25875453 0.32594526     0.051750906 0.06057729     0.28618251   0 0.13084112 0.04519758 0.2905217  0    0    0    1    3    0    0    3    0    0
##      BAEA BANS BAOR BARS BCCH BEKI BGGN BHCO BHVI BLJA BOBO BRCR BRTH BTNW BWHA BWWA CANG CARW CEDW CERW CHSP CHSW CLSW COGR COHA COYE DEJU DOWO EABL EAKI EAME EAPH EATO EAWP EUST FISP GBHE GCFL GHOW GRCA GRHE HAWO HERG HETH HOFI HOSP HOWA HOWR INBU KILL LOWA MALL MODO NOCA NOFL NOPA NRWS OROR OVEN
## 1002    0    0    0    0    1    0    1    2    0    1    0    0    0    0    0    1    0    0    1    0    0    1    0    0    0    0    0    2    0    0    0    0    3    2    0    0    0    0    0    0    0    0    0    0    0    0    3    0    2    0    0    0    0    3    0    0    0    0    0
## 1003    0    0    2    0    2    0    1    2    0    2    0    0    0    0    0    0    0    0    1    0    0    0    0    2    0    0    0    1    0    0    0    0    0    3    0    0    0    1    0    0    0    0    0    0    0    0    0    1    0    0    0    0    0    5    0    0    0    0    0
## 1004    0    0    0    0    6    1    2    3    0    6    0    0    0    0    0    0    1    0    0    0    0    2    0    0    0    0    1    3    0    0    0    0    0    7    0    0    0    1    0    0    0    1    0    0    0    0    0    1    2    0    0    0    0    3    1    0    0    0    0
## 1005    0    0    1    0    1    0    0    2    0    3    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    2    0    2    0    1    0    0    1    3    0    0    0    1    0    3    0    0    0    0    0    0    1    0    0    0    0    0    0    3    0    0    0    0    0
## 1006    0    0    1    0    6    0    2    0    0    7    0    0    0    0    0    0    0    0    0    0    0    1    0    1    0    0    0    5    0    0    0    0    0    7    0    0    0    3    0    0    0    3    0    0    0    0    0    0    0    0    0    0    0    5    1    0    0    0    0
## 1008    0    0    0    0    1    0    1    0    0    2    0    0    0    0    0    0    0    0    1    0    0    0    0    2    0    0    0    3    0    0    0    0    0    3    1    0    0    1    0    1    0    0    0    0    1    0    0    3    0    0    0    0    0    2    1    0    0    0    0
##      PIWO PUFI PUMA RBGR RBGU RBNU RBWO REVI RHWO ROPI RSHA RTHA RTHU RWBL SAVS SCTA SOSP SWSP TEWA TRES TUTI TUVU VEER VIRA WAVI WBNU WEVI WIFL WITU WIWR WODU WOTH YBCU YEWA YTVI YTWA
## 1002    0    0    0    0    0    0    0    2    0    0    0    0    0    0    0    0    0    0    0    0    1    0    0    0    0    1    0    0    0    0    0    2    0    0    1    0
## 1003    0    0    0    0    0    0    1    5    0    0    0    0    0    0    0    2    0    0    0    0    2    0    1    0    0    5    0    0    0    0    0    4    0    0    0    0
## 1004    0    0    0    0    2    0    5    7    0    0    1    0    0    0    0    5    0    0    0    0    6    0    0    0    0    4    0    0    0    0    0    3    0    0    0    0
## 1005    1    0    0    1    0    0    3    1    0    0    0    0    1    3    0    1    2    3    0    0    3    0    0    0    0    3    0    1    0    0    0    2    0    0    0    0
## 1006    4    0    0    4    0    0    6    6    0    0    0    0    0    0    0    6    0    0    0    0    4    0    2    0    0    3    0    0    0    1    0    6    0    0    0    0
## 1008    0    0    0    0    0    0    2    2    2    0    0    1    0    0    0    0    0    0    0    0    3    0    0    0    0    2    0    0    0    0    0    0    0    0    2    0

This is a pretty unwieldly dataset at present - we have 164 rows (representing our BBS/PCAP plots) and 131 columns. Out of those 131 columns, 104 are unique species (designated with four-letter codes), so most of the first 27 covariates are essentially our “predictors” - in this case, our PCAP data.

Let’s do a bit more processing on this dataset to get it in shape for downstream analyses. First, we will separate out our bird community data.

birdcomm<-CMP_bbs_pcap_join_pcaponly_pcapsum[,27:(ncol(CMP_bbs_pcap_join_pcaponly_pcapsum)-1)]

birdcomm[is.na(birdcomm)] <- 0

Next, we separate out our PCAP data.

vegecomm<-CMP_bbs_pcap_join_pcaponly_pcapsum[,4:25]

Since we are using the PCAP dataset as predictors here, we should likely check for correlations between these covariates. I have selected a (semi-)arbitrary cut-off of 0.6 to identify covariates as correlated; we then drop any covariates that are likely redundant.

vegecomm_cor<-cor(vegecomm)

set.seed(12345)
corvar<-findCorrelation(vegecomm_cor,cutoff=0.6)

vegecomm_drop<-vegecomm[,-corvar]

Now, we have:

  1. a vegecomm dataset - which representing PCAP measurements for each of the 164 BBS/PCAP plots.
  2. a birdcomm dataset - representing presence/absence of each bird species for each survey in each plot.

Note that we are not explicitly considering differences in bird counts (or occupancy, or detectability) through time at this point - we will get to that later. For now, we just want to explore broad associations between these two datasets!

Note as of 6/13/2022 - we will also want to only work with only binary data.

Community analyses

When working with datasets like these, we first have to acknowledge that these datasets are multivariate - with many possible predictor covariates (for example, plot-level characteristics) and many possible response covariates (for example, presence/absence of species in a plot). Fortunately, community ecology as a field has developed a number of standard approaches for exploring multivariate datasets, of which the most common is (arguably!) ordination. In essence, ordination techniques try to shrink down highly multidimensional datasets into a more manageable number of dimensions - often just a handful (and for initial plotting, typically just two). There are many, many possible ordination tools available for exploring community ecological datasets, and we will explore a few.

unconstrained ordination

In this section, we will be exploring ways of exploring the bird dataset on its own - basically, without any predictors.

First, we must reclassify our bird community data onto a binary scale - ranging from 0 (no detection) to 1 (detection).

birdcomm_bin<-birdcomm

birdcomm_bin<- birdcomm_bin %>% mutate_if(is.numeric, ~1 * (. >= 1))

One of the most common approaches for unconstrained ordination is called Principal Components Analysis (PCA). I have illustrated two approaches for this below - using either the prcomp() function or vegan’s rda() function (which we will use a bit more below!).

PCA.res<-prcomp(birdcomm_bin, scale=T)

biplot(PCA.res)
title("PCA (prcomp version)")

Personally, I find rda() to be much easier to use, so I tend to use it preferentially. Here is a slightly more complicated plot compared to the above - this time with convex hulls and points colored by reservation.

reserv<-factor(toupper(CMP_bbs_pcap_join_pcaponly_pcapsum$reservatio))
colpal_reserv<-rainbow(length(levels(reserv)))

p_unconc<-rda(birdcomm_bin,scale=T)

PCAscores <- scores(p_unconc, display = "sites") %>% 
  as.data.frame() %>% 
  rownames_to_column("site") 

PCAvect <- scores(p_unconc, display = "species") %>% 
  as.data.frame()

reservplot<-ordiplot(p_unconc,scaling=3,type="n")
points(reservplot, "sites",pch=21, col=colpal_reserv[reserv])


ordihull(reservplot,groups = reserv, col=colpal_reserv, lwd=3)
title("PCA (vegan rda version)")
legend(x="right", legend=levels(reserv), col=colpal_reserv,pch = 1,cex=0.5) 

while a bit hard to interpret, plots within Rocky River (RR) seem to be the most variable - and it might be that we see higher diversity in that area. We will see as we proceed if that is true.

metric multidimensional scaling (Principal Coordinates Analysis)

You can think of this as a distance-based analog to PCA. In this case, since we are using Bray-Curtis distances on binary data (presence/absence), this is equivalent to Sorensen’s index.

non-metric multidimensional scaling

bird_nmds<-metaMDS(birdcomm_bin,distance = "bray",trymax = 400,k=3)
## Run 0 stress 0.1751474 
## Run 1 stress 0.1761345 
## Run 2 stress 0.1732776 
## ... New best solution
## ... Procrustes: rmse 0.03981098  max resid 0.1508202 
## Run 3 stress 0.1743699 
## Run 4 stress 0.1738508 
## Run 5 stress 0.1725355 
## ... New best solution
## ... Procrustes: rmse 0.01742725  max resid 0.1440867 
## Run 6 stress 0.1733826 
## Run 7 stress 0.1753801 
## Run 8 stress 0.1727212 
## ... Procrustes: rmse 0.02129917  max resid 0.1476814 
## Run 9 stress 0.172957 
## ... Procrustes: rmse 0.02323107  max resid 0.1467514 
## Run 10 stress 0.1724776 
## ... New best solution
## ... Procrustes: rmse 0.00256073  max resid 0.01515115 
## Run 11 stress 0.1725733 
## ... Procrustes: rmse 0.00451041  max resid 0.03826016 
## Run 12 stress 0.172205 
## ... New best solution
## ... Procrustes: rmse 0.01955409  max resid 0.1465464 
## Run 13 stress 0.1732051 
## Run 14 stress 0.175053 
## Run 15 stress 0.1752042 
## Run 16 stress 0.1763912 
## Run 17 stress 0.1745713 
## Run 18 stress 0.1729667 
## Run 19 stress 0.1734391 
## Run 20 stress 0.1729283 
## Run 21 stress 0.1732458 
## Run 22 stress 0.1752081 
## Run 23 stress 0.1754885 
## Run 24 stress 0.1738273 
## Run 25 stress 0.1721794 
## ... New best solution
## ... Procrustes: rmse 0.003391018  max resid 0.03187639 
## Run 26 stress 0.173727 
## Run 27 stress 0.1755421 
## Run 28 stress 0.172803 
## Run 29 stress 0.1728229 
## Run 30 stress 0.1725033 
## ... Procrustes: rmse 0.02094271  max resid 0.1538802 
## Run 31 stress 0.173439 
## Run 32 stress 0.1751567 
## Run 33 stress 0.1724449 
## ... Procrustes: rmse 0.01094501  max resid 0.05449719 
## Run 34 stress 0.1722481 
## ... Procrustes: rmse 0.0153566  max resid 0.1417752 
## Run 35 stress 0.1727234 
## Run 36 stress 0.1748387 
## Run 37 stress 0.1735728 
## Run 38 stress 0.1752443 
## Run 39 stress 0.1733246 
## Run 40 stress 0.1726289 
## ... Procrustes: rmse 0.02159124  max resid 0.1466149 
## Run 41 stress 0.1729111 
## Run 42 stress 0.173677 
## Run 43 stress 0.1747492 
## Run 44 stress 0.1721845 
## ... Procrustes: rmse 0.004780107  max resid 0.04519247 
## Run 45 stress 0.1744348 
## Run 46 stress 0.1728829 
## Run 47 stress 0.1745146 
## Run 48 stress 0.1725974 
## ... Procrustes: rmse 0.008605116  max resid 0.05020094 
## Run 49 stress 0.1767774 
## Run 50 stress 0.1738603 
## Run 51 stress 0.1725219 
## ... Procrustes: rmse 0.01557841  max resid 0.1479282 
## Run 52 stress 0.173304 
## Run 53 stress 0.1724993 
## ... Procrustes: rmse 0.01543618  max resid 0.145884 
## Run 54 stress 0.1764419 
## Run 55 stress 0.1734277 
## Run 56 stress 0.1731943 
## Run 57 stress 0.1725213 
## ... Procrustes: rmse 0.01543772  max resid 0.1475148 
## Run 58 stress 0.1784585 
## Run 59 stress 0.1740317 
## Run 60 stress 0.172172 
## ... New best solution
## ... Procrustes: rmse 0.002653639  max resid 0.02017705 
## Run 61 stress 0.1759132 
## Run 62 stress 0.1729089 
## Run 63 stress 0.173989 
## Run 64 stress 0.1758112 
## Run 65 stress 0.1736825 
## Run 66 stress 0.1733895 
## Run 67 stress 0.17459 
## Run 68 stress 0.1728892 
## Run 69 stress 0.1736083 
## Run 70 stress 0.1750135 
## Run 71 stress 0.1750304 
## Run 72 stress 0.1725164 
## ... Procrustes: rmse 0.01549591  max resid 0.1468597 
## Run 73 stress 0.1737165 
## Run 74 stress 0.1742695 
## Run 75 stress 0.1736979 
## Run 76 stress 0.1724548 
## ... Procrustes: rmse 0.01669034  max resid 0.1456585 
## Run 77 stress 0.1728158 
## Run 78 stress 0.1743856 
## Run 79 stress 0.1754929 
## Run 80 stress 0.1726014 
## ... Procrustes: rmse 0.01458213  max resid 0.1469235 
## Run 81 stress 0.1729228 
## Run 82 stress 0.1724446 
## ... Procrustes: rmse 0.01699636  max resid 0.1432058 
## Run 83 stress 0.1738788 
## Run 84 stress 0.1733452 
## Run 85 stress 0.1733887 
## Run 86 stress 0.1741618 
## Run 87 stress 0.1749985 
## Run 88 stress 0.1753223 
## Run 89 stress 0.1772594 
## Run 90 stress 0.1743855 
## Run 91 stress 0.1755166 
## Run 92 stress 0.1759072 
## Run 93 stress 0.1752903 
## Run 94 stress 0.1725475 
## ... Procrustes: rmse 0.01843695  max resid 0.1509723 
## Run 95 stress 0.1730089 
## Run 96 stress 0.1743684 
## Run 97 stress 0.174357 
## Run 98 stress 0.1722625 
## ... Procrustes: rmse 0.01545736  max resid 0.1381 
## Run 99 stress 0.172896 
## Run 100 stress 0.1732551 
## Run 101 stress 0.1757211 
## Run 102 stress 0.1750701 
## Run 103 stress 0.1728711 
## Run 104 stress 0.1759084 
## Run 105 stress 0.1727175 
## Run 106 stress 0.1756873 
## Run 107 stress 0.1738645 
## Run 108 stress 0.1742867 
## Run 109 stress 0.1740164 
## Run 110 stress 0.1752175 
## Run 111 stress 0.1723468 
## ... Procrustes: rmse 0.01984014  max resid 0.1459213 
## Run 112 stress 0.1760383 
## Run 113 stress 0.1728909 
## Run 114 stress 0.1725117 
## ... Procrustes: rmse 0.01563108  max resid 0.1470691 
## Run 115 stress 0.1744335 
## Run 116 stress 0.1752269 
## Run 117 stress 0.1736964 
## Run 118 stress 0.175521 
## Run 119 stress 0.1727817 
## Run 120 stress 0.1727144 
## Run 121 stress 0.172672 
## ... Procrustes: rmse 0.02113729  max resid 0.1400382 
## Run 122 stress 0.1765034 
## Run 123 stress 0.1752396 
## Run 124 stress 0.1727298 
## Run 125 stress 0.1740477 
## Run 126 stress 0.172341 
## ... Procrustes: rmse 0.01964859  max resid 0.145326 
## Run 127 stress 0.174959 
## Run 128 stress 0.1732102 
## Run 129 stress 0.172554 
## ... Procrustes: rmse 0.01377573  max resid 0.1275104 
## Run 130 stress 0.1725232 
## ... Procrustes: rmse 0.01199473  max resid 0.1173779 
## Run 131 stress 0.1757938 
## Run 132 stress 0.1769059 
## Run 133 stress 0.1725009 
## ... Procrustes: rmse 0.01591585  max resid 0.1463755 
## Run 134 stress 0.1745496 
## Run 135 stress 0.172888 
## Run 136 stress 0.1762493 
## Run 137 stress 0.1735169 
## Run 138 stress 0.1725253 
## ... Procrustes: rmse 0.01301501  max resid 0.1177242 
## Run 139 stress 0.1755574 
## Run 140 stress 0.1748397 
## Run 141 stress 0.1741356 
## Run 142 stress 0.1765357 
## Run 143 stress 0.1731915 
## Run 144 stress 0.172731 
## Run 145 stress 0.1721813 
## ... Procrustes: rmse 0.002955712  max resid 0.02239333 
## Run 146 stress 0.1723639 
## ... Procrustes: rmse 0.02035249  max resid 0.1485191 
## Run 147 stress 0.1721948 
## ... Procrustes: rmse 0.005701654  max resid 0.04621034 
## Run 148 stress 0.1734274 
## Run 149 stress 0.1728666 
## Run 150 stress 0.1730735 
## Run 151 stress 0.1739775 
## Run 152 stress 0.1731944 
## Run 153 stress 0.1722683 
## ... Procrustes: rmse 0.01543467  max resid 0.1331563 
## Run 154 stress 0.173358 
## Run 155 stress 0.1736016 
## Run 156 stress 0.1761904 
## Run 157 stress 0.1725567 
## ... Procrustes: rmse 0.01479097  max resid 0.1455197 
## Run 158 stress 0.1746467 
## Run 159 stress 0.1744006 
## Run 160 stress 0.1742816 
## Run 161 stress 0.1724327 
## ... Procrustes: rmse 0.01706571  max resid 0.1450522 
## Run 162 stress 0.1748837 
## Run 163 stress 0.1743231 
## Run 164 stress 0.1736804 
## Run 165 stress 0.1753484 
## Run 166 stress 0.1738681 
## Run 167 stress 0.1735758 
## Run 168 stress 0.1741785 
## Run 169 stress 0.1755309 
## Run 170 stress 0.1735338 
## Run 171 stress 0.1722336 
## ... Procrustes: rmse 0.01601707  max resid 0.1399591 
## Run 172 stress 0.1723744 
## ... Procrustes: rmse 0.02093664  max resid 0.1545463 
## Run 173 stress 0.1723332 
## ... Procrustes: rmse 0.0195882  max resid 0.14724 
## Run 174 stress 0.175532 
## Run 175 stress 0.1727458 
## Run 176 stress 0.1732635 
## Run 177 stress 0.1759102 
## Run 178 stress 0.1740541 
## Run 179 stress 0.1721916 
## ... Procrustes: rmse 0.005767637  max resid 0.04636342 
## Run 180 stress 0.1749051 
## Run 181 stress 0.173469 
## Run 182 stress 0.1726598 
## ... Procrustes: rmse 0.02291106  max resid 0.1514561 
## Run 183 stress 0.1724743 
## ... Procrustes: rmse 0.0203851  max resid 0.1472383 
## Run 184 stress 0.1752857 
## Run 185 stress 0.1739626 
## Run 186 stress 0.1738131 
## Run 187 stress 0.1721772 
## ... Procrustes: rmse 0.001878673  max resid 0.01768344 
## Run 188 stress 0.1745026 
## Run 189 stress 0.1726038 
## ... Procrustes: rmse 0.02171574  max resid 0.1486159 
## Run 190 stress 0.1742865 
## Run 191 stress 0.1722344 
## ... Procrustes: rmse 0.01624271  max resid 0.1427985 
## Run 192 stress 0.1743934 
## Run 193 stress 0.1725121 
## ... Procrustes: rmse 0.01557617  max resid 0.1470384 
## Run 194 stress 0.1740988 
## Run 195 stress 0.1729496 
## Run 196 stress 0.172682 
## Run 197 stress 0.1746104 
## Run 198 stress 0.1758245 
## Run 199 stress 0.1762314 
## Run 200 stress 0.1729143 
## Run 201 stress 0.1745287 
## Run 202 stress 0.1725161 
## ... Procrustes: rmse 0.01567508  max resid 0.1477249 
## Run 203 stress 0.1751867 
## Run 204 stress 0.1733611 
## Run 205 stress 0.172937 
## Run 206 stress 0.1740611 
## Run 207 stress 0.1722337 
## ... Procrustes: rmse 0.01782391  max resid 0.1466053 
## Run 208 stress 0.1745937 
## Run 209 stress 0.173274 
## Run 210 stress 0.1748937 
## Run 211 stress 0.1739471 
## Run 212 stress 0.1734346 
## Run 213 stress 0.1742989 
## Run 214 stress 0.1733951 
## Run 215 stress 0.1731381 
## Run 216 stress 0.1722191 
## ... Procrustes: rmse 0.006453694  max resid 0.04636308 
## Run 217 stress 0.1744902 
## Run 218 stress 0.1755554 
## Run 219 stress 0.1726841 
## Run 220 stress 0.1727686 
## Run 221 stress 0.1722391 
## ... Procrustes: rmse 0.01774966  max resid 0.147681 
## Run 222 stress 0.1767903 
## Run 223 stress 0.1752083 
## Run 224 stress 0.1726425 
## ... Procrustes: rmse 0.01806157  max resid 0.1415587 
## Run 225 stress 0.174496 
## Run 226 stress 0.1729073 
## Run 227 stress 0.1756222 
## Run 228 stress 0.174881 
## Run 229 stress 0.1736295 
## Run 230 stress 0.1749679 
## Run 231 stress 0.1725654 
## ... Procrustes: rmse 0.01652435  max resid 0.1563714 
## Run 232 stress 0.1749454 
## Run 233 stress 0.175765 
## Run 234 stress 0.1723714 
## ... Procrustes: rmse 0.01981196  max resid 0.144888 
## Run 235 stress 0.1751377 
## Run 236 stress 0.1739592 
## Run 237 stress 0.1723483 
## ... Procrustes: rmse 0.01780819  max resid 0.1386382 
## Run 238 stress 0.1725032 
## ... Procrustes: rmse 0.01570548  max resid 0.1460301 
## Run 239 stress 0.1739827 
## Run 240 stress 0.1763483 
## Run 241 stress 0.1730713 
## Run 242 stress 0.1733056 
## Run 243 stress 0.1727 
## Run 244 stress 0.175973 
## Run 245 stress 0.1725987 
## ... Procrustes: rmse 0.01640206  max resid 0.1588538 
## Run 246 stress 0.1730129 
## Run 247 stress 0.175144 
## Run 248 stress 0.1779461 
## Run 249 stress 0.1783918 
## Run 250 stress 0.1729347 
## Run 251 stress 0.1743098 
## Run 252 stress 0.1723768 
## ... Procrustes: rmse 0.02059765  max resid 0.1504712 
## Run 253 stress 0.1733464 
## Run 254 stress 0.1748556 
## Run 255 stress 0.1733471 
## Run 256 stress 0.1754993 
## Run 257 stress 0.1759528 
## Run 258 stress 0.1753172 
## Run 259 stress 0.1747468 
## Run 260 stress 0.176772 
## Run 261 stress 0.1746528 
## Run 262 stress 0.173722 
## Run 263 stress 0.1721774 
## ... Procrustes: rmse 0.002415078  max resid 0.01975738 
## Run 264 stress 0.1724053 
## ... Procrustes: rmse 0.02043148  max resid 0.1515156 
## Run 265 stress 0.1738366 
## Run 266 stress 0.1754295 
## Run 267 stress 0.1729009 
## Run 268 stress 0.1729242 
## Run 269 stress 0.1775915 
## Run 270 stress 0.1727292 
## Run 271 stress 0.1728591 
## Run 272 stress 0.1759786 
## Run 273 stress 0.1724741 
## ... Procrustes: rmse 0.02040129  max resid 0.1473668 
## Run 274 stress 0.1733485 
## Run 275 stress 0.1736517 
## Run 276 stress 0.1741278 
## Run 277 stress 0.1721775 
## ... Procrustes: rmse 0.002352542  max resid 0.02004201 
## Run 278 stress 0.1742175 
## Run 279 stress 0.1725201 
## ... Procrustes: rmse 0.01383473  max resid 0.1180847 
## Run 280 stress 0.1727011 
## Run 281 stress 0.1725257 
## ... Procrustes: rmse 0.020634  max resid 0.1548589 
## Run 282 stress 0.174761 
## Run 283 stress 0.1754632 
## Run 284 stress 0.1751458 
## Run 285 stress 0.1724602 
## ... Procrustes: rmse 0.02056143  max resid 0.1475234 
## Run 286 stress 0.1734248 
## Run 287 stress 0.1726185 
## ... Procrustes: rmse 0.02165034  max resid 0.141876 
## Run 288 stress 0.1728716 
## Run 289 stress 0.1745318 
## Run 290 stress 0.1737329 
## Run 291 stress 0.1753969 
## Run 292 stress 0.1760232 
## Run 293 stress 0.1723744 
## ... Procrustes: rmse 0.0056035  max resid 0.03506882 
## Run 294 stress 0.1731454 
## Run 295 stress 0.1721988 
## ... Procrustes: rmse 0.004916334  max resid 0.04590163 
## Run 296 stress 0.173888 
## Run 297 stress 0.1763478 
## Run 298 stress 0.1722104 
## ... Procrustes: rmse 0.004017098  max resid 0.0320913 
## Run 299 stress 0.1725101 
## ... Procrustes: rmse 0.008473902  max resid 0.05683405 
## Run 300 stress 0.1755023 
## Run 301 stress 0.1726061 
## ... Procrustes: rmse 0.01799632  max resid 0.1455019 
## Run 302 stress 0.1727611 
## Run 303 stress 0.1736767 
## Run 304 stress 0.1728906 
## Run 305 stress 0.1740513 
## Run 306 stress 0.1735696 
## Run 307 stress 0.1739728 
## Run 308 stress 0.1787724 
## Run 309 stress 0.1736607 
## Run 310 stress 0.1737475 
## Run 311 stress 0.1727552 
## Run 312 stress 0.1726611 
## ... Procrustes: rmse 0.02117768  max resid 0.1404394 
## Run 313 stress 0.1733226 
## Run 314 stress 0.1727581 
## Run 315 stress 0.1732464 
## Run 316 stress 0.1723731 
## ... Procrustes: rmse 0.02080382  max resid 0.151506 
## Run 317 stress 0.1752211 
## Run 318 stress 0.175093 
## Run 319 stress 0.1766787 
## Run 320 stress 0.1721685 
## ... New best solution
## ... Procrustes: rmse 0.0009834757  max resid 0.01024635 
## Run 321 stress 0.1751463 
## Run 322 stress 0.1741857 
## Run 323 stress 0.1747149 
## Run 324 stress 0.1739038 
## Run 325 stress 0.1733386 
## Run 326 stress 0.1722203 
## ... Procrustes: rmse 0.004457349  max resid 0.03053609 
## Run 327 stress 0.175816 
## Run 328 stress 0.175337 
## Run 329 stress 0.1763445 
## Run 330 stress 0.1726706 
## Run 331 stress 0.1754963 
## Run 332 stress 0.1747553 
## Run 333 stress 0.173047 
## Run 334 stress 0.1739535 
## Run 335 stress 0.1747817 
## Run 336 stress 0.1724755 
## ... Procrustes: rmse 0.02017503  max resid 0.1470259 
## Run 337 stress 0.1758235 
## Run 338 stress 0.1735397 
## Run 339 stress 0.1721884 
## ... Procrustes: rmse 0.004215789  max resid 0.04444782 
## Run 340 stress 0.1738999 
## Run 341 stress 0.1769098 
## Run 342 stress 0.1751327 
## Run 343 stress 0.1730756 
## Run 344 stress 0.1746818 
## Run 345 stress 0.1725896 
## ... Procrustes: rmse 0.02081706  max resid 0.1467707 
## Run 346 stress 0.1739574 
## Run 347 stress 0.1746009 
## Run 348 stress 0.1737595 
## Run 349 stress 0.1745198 
## Run 350 stress 0.1730399 
## Run 351 stress 0.1733815 
## Run 352 stress 0.1735082 
## Run 353 stress 0.1748932 
## Run 354 stress 0.175114 
## Run 355 stress 0.1733985 
## Run 356 stress 0.1726582 
## ... Procrustes: rmse 0.01868993  max resid 0.1450401 
## Run 357 stress 0.1730416 
## Run 358 stress 0.1755821 
## Run 359 stress 0.176917 
## Run 360 stress 0.1772515 
## Run 361 stress 0.1732641 
## Run 362 stress 0.1726746 
## Run 363 stress 0.1750414 
## Run 364 stress 0.1745283 
## Run 365 stress 0.1731789 
## Run 366 stress 0.1770884 
## Run 367 stress 0.1742736 
## Run 368 stress 0.1757257 
## Run 369 stress 0.1749761 
## Run 370 stress 0.1729735 
## Run 371 stress 0.173791 
## Run 372 stress 0.1723384 
## ... Procrustes: rmse 0.0194334  max resid 0.1454182 
## Run 373 stress 0.1727315 
## Run 374 stress 0.1726655 
## ... Procrustes: rmse 0.02124793  max resid 0.1450788 
## Run 375 stress 0.1733572 
## Run 376 stress 0.1723717 
## ... Procrustes: rmse 0.01966895  max resid 0.144545 
## Run 377 stress 0.172977 
## Run 378 stress 0.1722662 
## ... Procrustes: rmse 0.003896939  max resid 0.02814835 
## Run 379 stress 0.1727404 
## Run 380 stress 0.1731728 
## Run 381 stress 0.1729195 
## Run 382 stress 0.1761529 
## Run 383 stress 0.1726321 
## ... Procrustes: rmse 0.02534648  max resid 0.1557111 
## Run 384 stress 0.174935 
## Run 385 stress 0.1727799 
## Run 386 stress 0.1742647 
## Run 387 stress 0.1742623 
## Run 388 stress 0.1723457 
## ... Procrustes: rmse 0.01991658  max resid 0.1496722 
## Run 389 stress 0.1722464 
## ... Procrustes: rmse 0.01773663  max resid 0.1480927 
## Run 390 stress 0.1761951 
## Run 391 stress 0.1759265 
## Run 392 stress 0.1750947 
## Run 393 stress 0.1727965 
## Run 394 stress 0.176099 
## Run 395 stress 0.1733906 
## Run 396 stress 0.1728845 
## Run 397 stress 0.1725288 
## ... Procrustes: rmse 0.0127141  max resid 0.1179979 
## Run 398 stress 0.1752787 
## Run 399 stress 0.1762001 
## Run 400 stress 0.1750601 
## *** No convergence -- monoMDS stopping criteria:
##    266: no. of iterations >= maxit
##    134: stress ratio > sratmax
plot(bird_nmds)

stressplot(bird_nmds)

bird species most associated with community differences

We can also look at which birds are most differentiated across communities for both the PCOA and NMDS ordinations. The below code essentially a) relates individual bird presence/absence to the final PCOA ordination, b) extracts p-values for associations with MDS1 and 2 (the first two ordination axes), and c) selects only those birds that are most differentiated at p less than or equal to 0.001. We will do so for the PCA, PCOA, and NMDS.

For all of these analyses, keep in mind that these are differentiated across plots, rather than reservations. We will explore a way to identify reservation-level differentiated species a litle bit below.

Here is the PCA:

fitpca <- envfit(p_unconc, birdcomm_bin, perm = 999) 

# extract p-values for each species
fitpca_pvals <- fitpca$vectors$pvals %>% 
  as.data.frame() %>% 
  rownames_to_column("species") %>% 
  dplyr::rename("pvals" = ".")

# extract coordinates for species, only keep species with p-val = 0.001
fitpca_spp <- fitpca %>% 
  scores(., display = "vectors") %>% 
  as.data.frame() %>% 
  rownames_to_column("species") %>% 
  full_join(., fitpca_pvals, by = "species") %>% 
  filter(pvals == 0.001)

47 species were identified as significantly contributing to variation between plots. These are my best attempts to break these down by typical migratory patterns for each species in Ohio (courtesy of www.allaboutbirds.org).

Migratory, breeding range in Ohio: Acadian Fly-Catcher, Baltimore Oriole, Barn Swallow, Blue-Gray Gnatcatcher, Brown Thrasher, Common Yellowthroat, Eastern Bluebird, Eastern Kingbird, Eastern Towhee, Eastern Wood-Pewee, Gray Catbird, Great-Crested Flycatcher, Hooded Warbler, House Wren, Indigo Bunting, Killdeer, Louisiana Water-Thrush, Orchard Oriole, Rose-breasted Grosbeak, Red-eyed Vireo, Scarlet Tanager, Tree Swallow, Warbling Vireo, Willow Flycatcher, Wood Thrush, Yellow Warbler, Yellow-Throated Vireo = 27 species

Year-round Ohio resident: American Crow, American Kestrel, Black-Capped Chickadee, Brown-Headed Cowbird, Cedar Waxwing, Common Grackle, Eastern Meadowlark, European Starling, Hairy Woodpecker, House Finch, House Sparrow, Mourning Dove, Northern Flicker, Pileated Woodpecker, Red-winged Blackbird, Song Sparrow, Tufted Titmouse, Turkey Vulture, White-breasted Nuthatch = 19 species

Migratory, non-breeding range in Ohio: Dark-eyed Junco = 1 species

We can repeat this procedure for the PCOA results too. Whereas PCA (above) was based on summarizing our bird community in the fewest number of axes that explain the most variance, PCOA is trying to maximize distance between plots. This is somewhat semantic in many cases, but will lead to differing representations of which birds are most associated with ordination axes.

Anyway! on to the PCOA:

# envfit() takes the output of metaMDS() and the species matrix you created
fitpcoa <- envfit(bird_pcoa, birdcomm_bin, perm = 999) 

# extract p-values for each species
fitpcoa_pvals <- fitpcoa$vectors$pvals %>% 
  as.data.frame() %>% 
  rownames_to_column("species") %>% 
  dplyr::rename("pvals" = ".")

# extract coordinates for species, only keep species with p-val = 0.001
fit_pcoa_spp <- fitpcoa %>% 
  scores(., display = "vectors") %>% 
  as.data.frame() %>% 
  rownames_to_column("species") %>% 
  full_join(., fitpcoa_pvals, by = "species") %>% 
  filter(pvals == 0.001)

So, now we have 40 species that contribute most to community distance between plots.

Migratory, breeding range in Ohio: Acadian Fly-Catcher, Baltimore Oriole, Barn Swallow, Blue-Gray Gnatcatcher, Brown Thrasher, Chimney Swift, Common Yellowthroat, Eastern Kingbird, Eastern Wood-Pewee, Gray Catbird, Hooded Warbler, House Wren, Indigo Bunting, Louisiana Water-Thrush, Orchard Oriole, Rose-breasted Grosbeak, Red-eyed Vireo, Scarlet Tanager, Tree Swallow, Veery, Warbling Vireo, Willow Flycatcher, Wood Thrush, Yellow Warbler, Yellow-Throated Vireo = 25 species

Year-round Ohio resident: American Robin, Brown-Headed Cowbird, Cedar Waxwing, Common Grackle, European Starling, Field Sparrow, House Finch, House Sparrow, Mourning Dove, Northern Flicker, Pileated Woodpecker, Red-winged Blackbird, Song Sparrow, White-breasted Nuthatch = 14 species

Migratory, non-breeding range in Ohio: Dark-eyed Junco = 1 species

and now the NMDS:

##    species      NMDS1       NMDS2 pvals
## 1     ACFL -0.6615159 -0.33993204 0.001
## 2     AMCR -0.1824643 -0.21687480 0.001
## 3     AMRO  0.1652093  0.29759255 0.001
## 4     BAOR  0.4715001 -0.14130150 0.001
## 5     BARS  0.3847375 -0.18007515 0.001
## 6     BCCH -0.1821493 -0.27828156 0.001
## 7     BGGN -0.1695426 -0.29758264 0.001
## 8     BHCO  0.2126973 -0.30699109 0.001
## 9     BRTH  0.2836090 -0.25445830 0.001
## 10    CEDW  0.4963877 -0.14194614 0.001
## 11    COGR  0.4168400  0.19899144 0.001
## 12    COYE  0.4305116 -0.45656653 0.001
## 13    DEJU -0.2774805 -0.27502180 0.001
## 14    EAKI  0.3000369 -0.15315332 0.001
## 15    EAWP -0.5238590 -0.16395515 0.001
## 16    EUST  0.5560226  0.08649803 0.001
## 17    FISP  0.1808023 -0.31446546 0.001
## 18    GRCA  0.5924405  0.01215275 0.001
## 19    HOFI  0.2812091  0.22451853 0.001
## 20    HOSP  0.2934945  0.26677074 0.001
## 21    HOWA -0.3074823 -0.54894517 0.001
## 22    HOWR  0.3033252  0.22647663 0.001
## 23    INBU  0.2723384 -0.47973161 0.001
## 24    LOWA -0.1413892 -0.35488881 0.001
## 25    MODO  0.4502863 -0.27715337 0.001
## 26    NOFL  0.3207516 -0.04470348 0.001
## 27    OROR  0.6432191 -0.04993506 0.001
## 28    PIWO -0.1199985 -0.42279736 0.001
## 29    RBGR  0.2600685 -0.34687602 0.001
## 30    REVI -0.4280763 -0.13293623 0.001
## 31    RWBL  0.6453788 -0.06557360 0.001
## 32    SCTA -0.4184239 -0.49848494 0.001
## 33    SOSP  0.6433413 -0.09766379 0.001
## 34    TRES  0.4269827 -0.11404675 0.001
## 35    WAVI  0.5914979  0.02827896 0.001
## 36    WBNU -0.3438214 -0.24357260 0.001
## 37    WIFL  0.3800642 -0.27789620 0.001
## 38    WOTH -0.1567631 -0.52087264 0.001
## 39    YEWA  0.6798322 -0.20922757 0.001
## 40    YTVI -0.1698554 -0.34639158 0.001

So, now we have 41 species that contribute most to community distance between plots.

Migratory, breeding range in Ohio: Acadian Fly-Catcher, Baltimore Oriole, Barn Swallow, Blue-Gray Gnatcatcher, Brown Thrasher, Common Yellowthroat, Eastern Kingbird, Eastern Wood-Pewee, Gray Catbird, Hooded Warbler, House Wren, Indigo Bunting, Louisiana Water-Thrush, Orchard Oriole, Rose-breasted Grosbeak, Red-eyed Vireo, Scarlet Tanager, Tree Swallow, Veery, Warbling Vireo, Willow Flycatcher, Wood Thrush, Yellow Warbler, Yellow-Throated Vireo

Year-round Ohio resident: American Robin, Black-Capped Chickadee, Brown-Headed Cowbird, Canada Goose, Cedar Waxwing, Common Grackle, European Starling, Field Sparrow, House Finch, House Sparrow, Mourning Dove, Northern Flicker, Pileated Woodpecker, Red-winged Blackbird, Song Sparrow, White-breasted Nuthatch

Migratory, non-breeding range in Ohio: Dark-eyed Junco

While these three sets of species are pretty similar, we can also take a look at which species are different across approaches.

outersect(fit_spp$species,fit_pcoa_spp$species,fitpca_spp$species)
##  [1] "CHSW" "EATO" "VEER" "AMGO" "AMKE" "CANG" "EABL" "EAME" "GCFL" "KILL" "TUVU"

So American Crow, American Kestrel, Canada Goose, Chimney Swift, Eastern Bluebird, Eastern Meadowlark, Easter Towhee, Great-Crested Flycatcher, Hairy Woodpecker, Killdeer, Tufted Titmouse, and Turkey Vulture. At least to me, nothing particularly sticks out as a commonality across these species, but it may make just as much sense to simply pool these species across approaches:

sort(unique(c(fit_spp$species,fit_pcoa_spp$species,fitpca_spp$species)))
##  [1] "ACFL" "AMCR" "AMGO" "AMKE" "AMRO" "BAOR" "BARS" "BCCH" "BGGN" "BHCO" "BRTH" "CANG" "CEDW" "CHSW" "COGR" "COYE" "DEJU" "EABL" "EAKI" "EAME" "EATO" "EAWP" "EUST" "FISP" "GCFL" "GRCA" "HAWO" "HOFI" "HOSP" "HOWA" "HOWR" "INBU" "KILL" "LOWA" "MODO" "NOFL" "NRWS" "OROR" "PIWO" "RBGR" "REVI" "RWBL"
## [43] "SCTA" "SOSP" "TRES" "TUVU" "VEER" "WAVI" "WBNU" "WIFL" "WOTH" "YEWA" "YTVI"

This gives us a list of 52 bird species that seem to be particularly differentiated/variable across plots. That is precisely half of the total number of species in the overall dataset.

As we proceed, we will begin layering on our environmental datasets using three main approaches: constrained ordination, linear models of species richness, and multi-species occupancy models.

constrained

Whereas the ordination approaches above were focused on unconstrained ordination - in this case, exploring how species are arranged between sample sites - there is an entire suite of approaches that address how environmental gradients shape community assemblages among sample sites. These are called constrained ordination approaches. I am personally a huge proponent of redundancy analysis (RDA) as a flexible constrained ordination approach. Taken from the GUSTA ME tutorial on RDA:

RDA can also be considered a constrained version of principal components analysis (PCA), wherein canonical axes - built from linear combinations of response variables - must also be linear combinations of the explanatory variables (i.e. fitted by [multiple linear regression]). The RDA approach generates one ordination in the space defined by the matrix of response variables and another in the space defined by the matrix of explanatory variables.

So, in essence, the ordination plot produced will depict an optimal arrangement of environmental axes (RDA1, 2, etc.) that maximally explains variation in the species dataset.

We are going to use stepwise selection to explore which PCAP covariates most explain plot-level differences.

p<-rda(birdcomm_bin~.,vegecomm_drop,scale=T)
fit<-envfit(p,vegecomm_drop,scale=T)


p_0<-rda(birdcomm_bin~1,vegecomm_drop,scale=T)
p_1<-rda(birdcomm_bin~.,vegecomm_drop,scale=T)

set.seed(12345)
p_step<- ordistep(p_0, scope = formula(p_1))
## 
## Start: birdcomm_bin ~ 1 
## 
##                   Df    AIC      F Pr(>F)   
## + tol_rel          1 756.52 6.6214  0.005 **
## + svp              1 759.73 3.3565  0.005 **
## + sens_rel         1 760.09 2.9905  0.005 **
## + bryo             1 760.14 2.9425  0.005 **
## + shrub            1 760.70 2.3844  0.005 **
## + dicot            1 760.72 2.3572  0.005 **
## + sens_rel_notree  1 760.89 2.1925  0.005 **
## + cyper            1 761.26 1.8221  0.005 **
## + canopy           1 761.41 1.6719  0.005 **
## + subcanopy        1 761.12 1.9574  0.010 **
## + inv              1 760.60 2.4826  0.060 . 
## + hydro            1 761.74 1.3399  0.100 . 
## + small_tree       1 761.87 1.2124  0.120   
## + hydro_rel        1 762.00 1.0866  0.290   
## + ap               1 762.19 0.8948  0.500   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: birdcomm_bin ~ tol_rel 
## 
##           Df    AIC      F Pr(>F)   
## - tol_rel  1 761.09 6.6214  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   Df    AIC      F Pr(>F)   
## + dicot            1 756.28 2.2181  0.005 **
## + svp              1 756.37 2.1250  0.005 **
## + cyper            1 756.72 1.7772  0.005 **
## + shrub            1 757.08 1.4239  0.040 * 
## + hydro            1 757.11 1.3912  0.050 * 
## + canopy           1 757.22 1.2845  0.055 . 
## + bryo             1 757.02 1.4847  0.070 . 
## + inv              1 756.58 1.9146  0.080 . 
## + small_tree       1 757.29 1.2147  0.130   
## + subcanopy        1 757.41 1.0943  0.310   
## + hydro_rel        1 757.47 1.0400  0.345   
## + sens_rel_notree  1 757.46 1.0427  0.395   
## + sens_rel         1 757.50 1.0083  0.435   
## + ap               1 757.66 0.8463  0.600   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: birdcomm_bin ~ tol_rel + dicot 
## 
##           Df    AIC      F Pr(>F)   
## - dicot    1 756.52 2.2181  0.005 **
## - tol_rel  1 760.72 6.4527  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   Df    AIC      F Pr(>F)   
## + svp              1 756.32 1.9190  0.005 **
## + cyper            1 756.36 1.8840  0.005 **
## + bryo             1 756.74 1.5059  0.060 . 
## + inv              1 756.33 1.9124  0.075 . 
## + canopy           1 756.95 1.3002  0.080 . 
## + small_tree       1 757.04 1.2114  0.140   
## + shrub            1 757.17 1.0897  0.285   
## + sens_rel_notree  1 757.22 1.0397  0.340   
## + subcanopy        1 757.23 1.0279  0.390   
## + hydro_rel        1 757.25 1.0092  0.395   
## + sens_rel         1 757.23 1.0254  0.410   
## + ap               1 757.40 0.8590  0.550   
## + hydro            1 757.35 0.9127  0.695   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: birdcomm_bin ~ tol_rel + dicot + svp 
## 
##           Df    AIC      F Pr(>F)   
## - dicot    1 756.37 2.0114  0.010 **
## - svp      1 756.28 1.9190  0.005 **
## - tol_rel  1 759.74 5.3757  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   Df    AIC      F Pr(>F)   
## + cyper            1 756.50 1.7813  0.005 **
## + canopy           1 756.97 1.3150  0.070 . 
## + inv              1 756.35 1.9287  0.075 . 
## + bryo             1 756.98 1.3094  0.165   
## + small_tree       1 757.12 1.1738  0.200   
## + hydro_rel        1 757.27 1.0211  0.340   
## + shrub            1 757.29 1.0049  0.405   
## + subcanopy        1 757.29 1.0052  0.425   
## + sens_rel_notree  1 757.30 0.9966  0.440   
## + ap               1 757.43 0.8715  0.535   
## + sens_rel         1 757.34 0.9545  0.605   
## + hydro            1 757.38 0.9145  0.610   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: birdcomm_bin ~ tol_rel + dicot + svp + cyper 
## 
##           Df    AIC      F Pr(>F)   
## - cyper    1 756.32 1.7813  0.010 **
## - svp      1 756.36 1.8160  0.005 **
## - dicot    1 756.83 2.2826  0.005 **
## - tol_rel  1 759.99 5.4185  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   Df    AIC      F Pr(>F)  
## + inv              1 756.50 1.9396  0.080 .
## + bryo             1 757.12 1.3309  0.085 .
## + canopy           1 757.18 1.2729  0.100 .
## + small_tree       1 757.31 1.1473  0.260  
## + shrub            1 757.41 1.0509  0.320  
## + hydro_rel        1 757.41 1.0481  0.385  
## + sens_rel_notree  1 757.45 1.0108  0.400  
## + subcanopy        1 757.47 0.9963  0.525  
## + sens_rel         1 757.51 0.9578  0.560  
## + hydro            1 757.54 0.9208  0.600  
## + ap               1 757.70 0.7733  0.690  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_step
## Call: rda(formula = birdcomm_bin ~ tol_rel + dicot + svp + cyper, data = vegecomm_drop, scale = T)
## 
##                 Inertia Proportion Rank
## Total         103.00000    1.00000     
## Constrained     7.61479    0.07393    4
## Unconstrained  95.38521    0.92607  102
## Inertia is correlations 
## 
## Eigenvalues for constrained axes:
##  RDA1  RDA2  RDA3  RDA4 
## 4.483 1.406 0.991 0.735 
## 
## Eigenvalues for unconstrained axes:
##   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8 
## 6.458 3.725 3.535 3.198 3.030 2.685 2.578 2.387 
## (Showing 8 of 102 unconstrained eigenvalues)
anova(p_step)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = birdcomm_bin ~ tol_rel + dicot + svp + cyper, data = vegecomm_drop, scale = T)
##           Df Variance      F Pr(>F)    
## Model      4    7.615 3.1733  0.001 ***
## Residual 159   95.385                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(p_step,by="term")
## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = birdcomm_bin ~ tol_rel + dicot + svp + cyper, data = vegecomm_drop, scale = T)
##           Df Variance      F Pr(>F)    
## tol_rel    1    4.045 6.7420  0.001 ***
## dicot      1    1.345 2.2417  0.001 ***
## svp        1    1.157 1.9283  0.003 ** 
## cyper      1    1.069 1.7813  0.006 ** 
## Residual 159   95.385                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(p_step,by="axis")
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = birdcomm_bin ~ tol_rel + dicot + svp + cyper, data = vegecomm_drop, scale = T)
##           Df Variance      F Pr(>F)    
## RDA1       1    4.483 7.4724  0.001 ***
## RDA2       1    1.406 2.3436  0.001 ***
## RDA3       1    0.991 1.6527  0.007 ** 
## RDA4       1    0.735 1.2246  0.095 .  
## Residual 159   95.385                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ordiplot(p_step,scaling=3)

So this tells us that the most explanatory group of PCAP covariates is relative cover of tolerant plants (tol_rel), cover of dicots, cover of seedless vascular plants (svp), and cover of sedges (cyper). Rather conveniently, the relative cover of dicots, seedless vascular plants, and sedges is negatively associated with the relative cover of tolerant plant species. As a quick caveat, these four variables only explain about 7.3% of the variation, as seen here (red is our environmental axes, black is our unconstrained axes):

constrained_eig <- p_step$CCA$eig/p_step$tot.chi*100
unconstrained_eig <- p_step$CA$eig/p_step$tot.chi*100
expl_var <- c(constrained_eig, unconstrained_eig)
barplot (expl_var[1:20], col = c(rep ('red', length (constrained_eig)), rep ('black', length (unconstrained_eig))),
         las = 2, ylab = '% variation')

This is not that terribly uncommon for redundancy analyses - It is often difficult to find environmental gradients that explain all of the variance in large datasets like ours. Biologically, we likely do not necessarily expect that a single environmental gradient alone drives all of the variation in bird presence/absence, and intuitively know that we also have to worry about observer effects, weather, and other unmeasured characteristics. Regardless of any of the above, the ordination is significant - as are the relationships with these covariates, so we can safely say that nearly 10% of the variation in our bird species dataset is explainable by these 4 PCAP covariates!

We can also explore bird species most significantly associated with our constrained ordination axes, too.

# envfit() takes the output of metaMDS() and the species matrix you created
fitconcrda <- envfit(p_step, birdcomm_bin, perm = 999) 

# extract p-values for each species
fitconcrda_pvals <- fitconcrda$vectors$pvals %>% 
  as.data.frame() %>% 
  rownames_to_column("species") %>% 
  dplyr::rename("pvals" = ".")

# extract coordinates for species, only keep species with p-val = 0.001
fitconcrda_spp <- fitconcrda %>% 
  scores(., display = "vectors") %>% 
  as.data.frame() %>% 
  rownames_to_column("species") %>% 
  full_join(., fitconcrda_pvals, by = "species") %>% 
  filter(pvals == 0.001)

OK - this is getting towards our main analytical goals! Let’s sort by scores for that first environmental axis and plot these in order:

fitconcrda_spp_order<-fitconcrda_spp[order(fitconcrda_spp$RDA1),]
fitconcrda_spp_order$species<-factor(fitconcrda_spp_order$species,levels=fitconcrda_spp_order$species)


ggplot(aes(x=RDA1,y=species),data=fitconcrda_spp_order)+geom_point()+theme_classic(base_size=16)

This is a pretty neat plot - Acadian Flycatcher, Scarlet Tanager, Eastern Wood-Pewee, White-breasted Nuthatch, Hooded Warbler, Red-Eyed Vireo, Dark-Eyed Junco, Wood Thrush, Hairy Woodpeckers, and Lousiana Water-Thrush all appear to be associated with higher sedge, seedless vascular plant, and dicot cover and lower cover of tolerant plants. In contrast, Yellow Warblers, Orchard Orioles, Red-Winged Blackbirds, European Starlings, Warbling Vireos, Song Sparrows, Grey Catbirds, Tree Swallows, Mourning Doves, and Cedar Waxwings are all associated with more tolerant plants.

constrained - sites

While we have already explored environmental associations, we might also want to know what species are most differentiated between reservations. We will then use similar code to that used above to explore species that are most differentiated between reservations.

p_reservations<-rda(birdcomm_bin~reservatio,CMP_bbs_pcap_join_pcaponly_pcapsum,scale=T)

ordiplot(p_reservations,scaling=3)

reservations_fitbirds<-envfit(p_reservations,birdcomm_bin)

reservations_fit<-envfit(p_reservations~reservatio,data=CMP_bbs_pcap_join_pcaponly_pcapsum,scale=T)


CMP_bbs_pcap_join_pcaponly_pcapsum$site<-rownames(CMP_bbs_pcap_join_pcaponly_pcapsum)

plot_df_reservations <- scores(p_reservations, display = "sites") %>% 
  as.data.frame() %>% 
  rownames_to_column("site") %>% 
  full_join(CMP_bbs_pcap_join_pcaponly_pcapsum, by = "site")

plot_rda_reservations <- ggplot(plot_df_reservations, aes(x = RDA1, y = RDA2, color = reservatio)) +
  geom_point(size = 3, alpha = 0.8) +
  stat_ellipse(linetype = 2, size = 1) +
  labs(title = "RDA")
plot_rda_reservations
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Warning: Removed 7 row(s) containing missing values (geom_path).

fitreservations_pvals <- reservations_fitbirds$vectors$pvals %>% 
  as.data.frame() %>% 
  rownames_to_column("species") %>% 
  dplyr::rename("pvals" = ".")

# extract coordinates for species, only keep species with p-val = 0.001
fit_reservations_spp <- reservations_fitbirds %>% 
  scores(., display = "vectors") %>% 
  as.data.frame() %>% 
  rownames_to_column("species") %>% 
  full_join(., fitreservations_pvals, by = "species") %>% 
  filter(pvals == 0.001)
fit_reservations_spp
##    species        RDA1        RDA2 pvals
## 1     ACFL -0.67536570  0.09819212 0.001
## 2     AMCR -0.37692846 -0.20456090 0.001
## 3     AMRE -0.07400135 -0.35850489 0.001
## 4     AMRO  0.26495737  0.23369097 0.001
## 5     BAOR  0.41712060  0.17739155 0.001
## 6     BARS  0.36240761 -0.30012415 0.001
## 7     BCCH -0.24236822  0.20415669 0.001
## 8     BRCR -0.28606094 -0.09945096 0.001
## 9     BRTH  0.25037094 -0.24226047 0.001
## 10    CEDW  0.41740795 -0.28403651 0.001
## 11    COGR  0.51954277  0.18156226 0.001
## 12    COYE  0.09098344 -0.54427390 0.001
## 13    DEJU -0.37988349 -0.03903739 0.001
## 14    EATO -0.12707725 -0.27487048 0.001
## 15    EAWP -0.44389257  0.47561627 0.001
## 16    EUST  0.55101057 -0.40645368 0.001
## 17    FISP -0.02103263 -0.34546312 0.001
## 18    GRCA  0.46149939 -0.09909873 0.001
## 19    HAWO -0.14783078  0.26640131 0.001
## 20    HOFI  0.39273807 -0.02749710 0.001
## 21    HOSP  0.47906450  0.11992733 0.001
## 22    HOWA -0.56788748 -0.25735546 0.001
## 23    HOWR  0.32157965  0.13443633 0.001
## 24    INBU  0.14260564 -0.26880856 0.001
## 25    LOWA -0.29872249 -0.25058885 0.001
## 26    MODO  0.36998816 -0.15778881 0.001
## 27    NOFL  0.30988119  0.16027358 0.001
## 28    NRWS  0.44336152  0.10254562 0.001
## 29    OROR  0.54726150 -0.40264845 0.001
## 30    OVEN -0.27119395 -0.22211081 0.001
## 31    RBGU  0.36218342 -0.09026430 0.001
## 32    REVI -0.29178905  0.32096716 0.001
## 33    RWBL  0.50726734 -0.28022387 0.001
## 34    SCTA -0.54059278  0.06735648 0.001
## 35    SOSP  0.46731737 -0.16498066 0.001
## 36    TRES  0.44251225 -0.11105510 0.001
## 37    VEER -0.25041023 -0.17876454 0.001
## 38    WAVI  0.54355113 -0.29016743 0.001
## 39    WBNU -0.45277798  0.32182366 0.001
## 40    WIFL  0.28267826 -0.33559442 0.001
## 41    WOTH -0.35735577 -0.01891986 0.001
## 42    YEWA  0.50674786 -0.44480078 0.001
## 43    YTVI -0.35813458 -0.21812489 0.001
fit_reservations <- reservations_fit %>% 
  scores(., display = "factors") %>% 
  as.data.frame() %>%
  rownames_to_column("reservations")

rda_reservation_plot_new <- ggplot(plot_df_reservations, aes(x = RDA1, y = RDA2)) +
  coord_fixed() +
  geom_point(aes(color = reservatio), size = 3, alpha = 0.8) +
  geom_segment(data = fit_reservations_spp, aes(x = 0, xend = RDA1, y = 0, yend = RDA2),
               col = "black") +
  geom_text(data = fit_reservations_spp, aes(label = species))+
  geom_text(data=as.data.frame(fit_reservations),aes(label=reservations))
rda_reservation_plot_new

Interesting! So when we consider reservation-of-origin, Hinckley seems to have a pretty unique bird community overall, and groups negatively along RDA1. In contrast, Brookside, Ohio and Erie Canal, and Washington all group positively along RDA1 - and are all quite heavily urbanized.

While I personally think the environmental gradient analyses are a bit more interesting, this does highlight reservations that are more or less distinct in bird community space. We can also plot just the reservations alone too:

ggplot(plot_df_reservations, aes(x = RDA1, y = RDA2)) +
  coord_fixed() +
  geom_text(data=as.data.frame(fit_reservations),aes(label=reservations))+xlim(-2,4)+ylim(-4,2)

Now pRDA - with variance partitioning.

As a final sanity check, it can sometimes be difficult to disentangle geographic effects from environmental effects in these sorts of analyses. For example, what if (hypothetically) cover of seedless vascular plants is negatively associated with geographic coordinates? Thus, we might want to briefly check if the variance explained by our environmental axes is (or is not) confounded with geography. To do so, we will use the varpart() function, which partitions variance between model components:

vegecomm_latlong<-CMP_bbs_pcap_join_pcaponly_pcapsum[,2:26]

pcap_geog<-varpart(birdcomm_bin,~ tol_rel + svp + dicot + cyper,~longitude+latitude,data=vegecomm_latlong)
plot(pcap_geog)

So - most of this variation is associated with these top covariates, and NOT with latitude and longitude alone. Good to see - this means we can be confident that these are actual environmental effects!

community richness metrics

Next, we can explore the relationships between species richness and plot characteristics. We will start with plot-level explorations.

plot-level analyses

plot-level richness and various indicators

Now, we will produce some exploratory models and plots of richness against basic indicators of vegetative quality, including the Floristic Quality Indicator (FQI) and the Vegetative Index of Biotic Integrity (VIBI).

S <- specnumber(birdcomm_bin)


vegecomm_richness<-cbind(vegecomm,S)

cor(vegecomm_richness)
##                         vibifq vibifqnotrees        carex        cyper       dicot        shade       shrub        hydro         svp           ap           fqi         bryo   hydro_rel hydro_rel_notrees     sens_rel sens_rel_notree     tol_rel tol_rel_notree          inv    small_tree    subcanopy
## vibifq             1.000000000   0.903466506  0.302129838  0.268662618  0.44969374  0.643467211  0.42091385  0.094920118  0.52435545 -0.201584848  9.191368e-01  0.227249850  0.14141059       0.232333080  0.615943742     0.517291454 -0.82095096   -0.815660534 -0.204274958 -8.896552e-03  0.352987107
## vibifqnotrees      0.903466506   1.000000000  0.315594820  0.286141876  0.57241339  0.762505137  0.44919931  0.183868172  0.57125442 -0.163682237  9.157120e-01  0.179820785  0.23974361       0.250108233  0.401409023     0.573829627 -0.66469643   -0.758049585 -0.150119819  4.913246e-03  0.428733969
## carex              0.302129838   0.315594820  1.000000000  0.990155702  0.46569343  0.543313784  0.14791367  0.173971056  0.32038084 -0.207379567  4.247642e-01 -0.088332875 -0.01680755      -0.001503659 -0.010702821     0.022201528 -0.08217802   -0.040658648 -0.068732762  4.256074e-02  0.098208268
## cyper              0.268662618   0.286141876  0.990155702  1.000000000  0.45810906  0.519855636  0.12970354  0.173572197  0.30826804 -0.214455749  3.923608e-01 -0.092992405 -0.01088781       0.001649045 -0.026296645     0.010052307 -0.04666473   -0.008993867 -0.072806020  4.618045e-02  0.105018429
## dicot              0.449693738   0.572413388  0.465693427  0.458109060  1.00000000  0.900954560  0.38443818  0.501981435  0.48161327 -0.120998160  6.868442e-01 -0.078551692  0.22182199       0.148547049 -0.038371772     0.153124228 -0.09258616   -0.111832162 -0.074575448 -1.060284e-01  0.325216371
## shade              0.643467211   0.762505137  0.543313784  0.519855636  0.90095456  1.000000000  0.40882995  0.378278862  0.58861737 -0.156097164  8.368920e-01  0.020607143  0.20374398       0.157734262  0.086603981     0.266173043 -0.31108544   -0.344672034 -0.059325341  8.208519e-03  0.371162537
## shrub              0.420913850   0.449199309  0.147913674  0.129703536  0.38443818  0.408829955  1.00000000  0.588835377  0.29470052 -0.101698461  4.747080e-01  0.085207890  0.49360022       0.536538046  0.162917290     0.160089484 -0.22690227   -0.273227137 -0.107336552 -8.936686e-02  0.255636573
## hydro              0.094920118   0.183868172  0.173971056  0.173572197  0.50198144  0.378278862  0.58883538  1.000000000  0.17700645 -0.039802022  2.285981e-01 -0.081936068  0.54679684       0.488784912 -0.086929406    -0.061281690  0.10928036    0.079007210  0.003619817 -3.202940e-02  0.204315768
## svp                0.524355447   0.571254416  0.320380841  0.308268040  0.48161327  0.588617368  0.29470052  0.177006451  1.00000000 -0.147271791  5.999081e-01  0.244051454  0.15150342       0.153627835  0.190425176     0.258147075 -0.27430569   -0.331328168 -0.113257408 -9.341899e-02  0.268499329
## ap                -0.201584848  -0.163682237 -0.207379567 -0.214455749 -0.12099816 -0.156097164 -0.10169846 -0.039802022 -0.14727179  1.000000000 -2.573751e-01  0.006014263 -0.03413775      -0.021095091 -0.055016524    -0.012339374  0.12954362    0.125069560  0.033126004 -7.653394e-02 -0.095063187
## fqi                0.919136798   0.915712021  0.424764247  0.392360839  0.68684421  0.836892017  0.47470799  0.228598128  0.59990814 -0.257375101  1.000000e+00  0.164507181  0.15874283       0.185026315  0.395205020     0.452238507 -0.64079294   -0.658950801 -0.158986465 -1.659484e-05  0.356068079
## bryo               0.227249850   0.179820785 -0.088332875 -0.092992405 -0.07855169  0.020607143  0.08520789 -0.081936068  0.24405145  0.006014263  1.645072e-01  1.000000000  0.02837391       0.046954292  0.289335075     0.286497590 -0.32259935   -0.349022222 -0.101490062  3.167128e-01 -0.004418917
## hydro_rel          0.141410593   0.239743608 -0.016807547 -0.010887815  0.22182199  0.203743981  0.49360022  0.546796836  0.15150342 -0.034137749  1.587428e-01  0.028373905  1.00000000       0.899920418 -0.091789549    -0.064596088 -0.12854831   -0.192750547 -0.048956575  6.854047e-02  0.412103596
## hydro_rel_notrees  0.232333080   0.250108233 -0.001503659  0.001649045  0.14854705  0.157734262  0.53653805  0.488784912  0.15362783 -0.021095091  1.850263e-01  0.046954292  0.89992042       1.000000000  0.061436676    -0.073288844 -0.24903476   -0.272482935 -0.076359781  2.939615e-02  0.387893166
## sens_rel           0.615943742   0.401409023 -0.010702821 -0.026296645 -0.03837177  0.086603981  0.16291729 -0.086929406  0.19042518 -0.055016524  3.952050e-01  0.289335075 -0.09178955       0.061436676  1.000000000     0.529677884 -0.55933338   -0.537226545 -0.075091400 -7.745522e-03  0.003765294
## sens_rel_notree    0.517291454   0.573829627  0.022201528  0.010052307  0.15312423  0.266173043  0.16008948 -0.061281690  0.25814708 -0.012339374  4.522385e-01  0.286497590 -0.06459609      -0.073288844  0.529677884     1.000000000 -0.38059052   -0.445572344  0.002392817  2.451389e-02  0.060370376
## tol_rel           -0.820950963  -0.664696430 -0.082178015 -0.046664726 -0.09258616 -0.311085443 -0.22690227  0.109280360 -0.27430569  0.129543623 -6.407929e-01 -0.322599345 -0.12854831      -0.249034760 -0.559333375    -0.380590517  1.00000000    0.923531143  0.266063420 -1.012597e-01 -0.288508481
## tol_rel_notree    -0.815660534  -0.758049585 -0.040658648 -0.008993867 -0.11183216 -0.344672034 -0.27322714  0.079007210 -0.33132817  0.125069560 -6.589508e-01 -0.349022222 -0.19275055      -0.272482935 -0.537226545    -0.445572344  0.92353114    1.000000000  0.214269092 -9.522446e-02 -0.323032669
## inv               -0.204274958  -0.150119819 -0.068732762 -0.072806020 -0.07457545 -0.059325341 -0.10733655  0.003619817 -0.11325741  0.033126004 -1.589865e-01 -0.101490062 -0.04895657      -0.076359781 -0.075091400     0.002392817  0.26606342    0.214269092  1.000000000  6.041087e-02 -0.073180439
## small_tree        -0.008896552   0.004913246  0.042560741  0.046180446 -0.10602839  0.008208519 -0.08936686 -0.032029403 -0.09341899 -0.076533935 -1.659484e-05  0.316712780  0.06854047       0.029396146 -0.007745522     0.024513894 -0.10125967   -0.095224459  0.060410871  1.000000e+00 -0.150833647
## subcanopy          0.352987107   0.428733969  0.098208268  0.105018429  0.32521637  0.371162537  0.25563657  0.204315768  0.26849933 -0.095063187  3.560681e-01 -0.004418917  0.41210360       0.387893166  0.003765294     0.060370376 -0.28850848   -0.323032669 -0.073180439 -1.508336e-01  1.000000000
## canopy             0.203002590   0.102153956 -0.102731833 -0.111819839 -0.02406017  0.020867114  0.10411367  0.040288929  0.07700871 -0.055629052  1.434842e-01  0.104404368  0.04753125       0.107445181  0.264437686     0.096605726 -0.30995440   -0.253543359 -0.102465127  5.548526e-02  0.015482845
## S                 -0.180202756  -0.094422012 -0.123628799 -0.116601889  0.14178910  0.010138092 -0.09358624  0.123776993 -0.11846054  0.087142662 -1.110069e-01 -0.190998403 -0.03126316      -0.103537313 -0.096502774    -0.027700516  0.23555710    0.198107496  0.047266522  2.431789e-02 -0.094543020
##                        canopy           S
## vibifq             0.20300259 -0.18020276
## vibifqnotrees      0.10215396 -0.09442201
## carex             -0.10273183 -0.12362880
## cyper             -0.11181984 -0.11660189
## dicot             -0.02406017  0.14178910
## shade              0.02086711  0.01013809
## shrub              0.10411367 -0.09358624
## hydro              0.04028893  0.12377699
## svp                0.07700871 -0.11846054
## ap                -0.05562905  0.08714266
## fqi                0.14348423 -0.11100691
## bryo               0.10440437 -0.19099840
## hydro_rel          0.04753125 -0.03126316
## hydro_rel_notrees  0.10744518 -0.10353731
## sens_rel           0.26443769 -0.09650277
## sens_rel_notree    0.09660573 -0.02770052
## tol_rel           -0.30995440  0.23555710
## tol_rel_notree    -0.25354336  0.19810750
## inv               -0.10246513  0.04726652
## small_tree         0.05548526  0.02431789
## subcanopy          0.01548284 -0.09454302
## canopy             1.00000000 -0.08082640
## S                 -0.08082640  1.00000000
S_vibi<-lm(S~vibifq,vegecomm_richness)
summary(S_vibi)
## 
## Call:
## lm(formula = S ~ vibifq, data = vegecomm_richness)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.0444  -4.2186  -0.2331   3.8289  19.0802 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 28.41399    1.07526  26.425   <2e-16 ***
## vibifq      -0.04908    0.02105  -2.332   0.0209 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.217 on 162 degrees of freedom
## Multiple R-squared:  0.03247,    Adjusted R-squared:  0.0265 
## F-statistic: 5.437 on 1 and 162 DF,  p-value: 0.02094
S_fqi<-lm(S~fqi,vegecomm_richness)
summary(S_fqi)
## 
## Call:
## lm(formula = S ~ fqi, data = vegecomm_richness)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.1875  -4.4209  -0.1052   3.9411  19.4469 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.1082     1.4444  19.461   <2e-16 ***
## fqi          -0.1079     0.0759  -1.422    0.157    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.282 on 162 degrees of freedom
## Multiple R-squared:  0.01232,    Adjusted R-squared:  0.006226 
## F-statistic: 2.021 on 1 and 162 DF,  p-value: 0.157
ggplot(aes(x=vibifq,y=S),data=vegecomm_richness)+geom_point()+geom_smooth()

ggplot(aes(x=fqi,y=S),data=vegecomm_richness)+geom_point()+geom_smooth()

Somewhat unexpectedly, this seems to suggest that higher VIBI scores actually produce slightly less species-rich communities. What could be causing this? Well, to start - perhaps the number of surveys for each plot may be influencing this.

We can generate a summary of the number of surveys for each plot.

pcap_surveynumber<-CMP_bbs_pcap_join_pcaponly %>% 
  group_by(simp_plot) %>%
  summarise(n=n_distinct(date)
            ) 
vegecomm_richness$surveycount<-pcap_surveynumber$n

specaccum(birdcomm_bin,w=pcap_surveynumber$n,method="random")
## Species Accumulation Curve
## Accumulation method: random, with 100 permutations, weighted
## Call: specaccum(comm = birdcomm_bin, method = "random", w = pcap_surveynumber$n) 
## 
##                                                                                                                                                                                                                                                                                                    
## Sites    1.0000  2.0000  3.0000  4.0000  5.0000  6.0000  7.0000  8.0000  9.0000 10.0000 11.000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000 21.0000  22.0000  23.0000  24.0000  25.0000  26.0000  27.0000  28.0000  29.0000  30.0000  31.0000  32.0000  33.0000  34.0000
## Effort   4.6646  9.3293 13.9939 18.6585 23.3232 27.9878 32.6524 37.3171 41.9817 46.6463 51.311 55.9756 60.6402 65.3049 69.9695 74.6341 79.2988 83.9634 88.6280 93.2927 97.9573 102.6220 107.2866 111.9512 116.6159 121.2805 125.9451 130.6098 135.2744 139.9390 144.6037 149.2683 153.9329 158.5976
## Richness     NA 37.0601 43.9160 48.4527 52.0466 55.2858 58.1210 60.0900 61.7656 63.1583 64.629 66.1469 67.4490 68.7554 69.9239 70.9281 71.7678 72.6455 73.4860 74.3025 74.8494  75.6023  76.3352  76.8898  77.5041  77.9124  78.4497  79.0296  79.4300  79.9113  80.3283  80.6705  81.0429  81.4608
## sd           NA  5.4585  5.3625  5.5323  5.7458  5.6377  5.4931  5.4603  5.1609  5.0987  4.744  4.5960  4.5160  4.2286  4.2348  4.1589  4.0009  4.1083  4.1018  3.9875  4.0215   3.9500   3.7498   3.6637   3.6386   3.5648   3.5230   3.3841   3.3204   3.2723   3.2496   3.2391   3.1904   3.1615
##                                                                                                                                                                                                                                                                                                         
## Sites     35.0000  36.0000  37.0000  38.0000  39.0000  40.0000  41.0000  42.0000  43.0000  44.0000  45.0000  46.0000  47.0000  48.0000  49.0000  50.0000  51.0000  52.0000  53.0000  54.0000  55.0000  56.0000  57.0000  58.0000  59.0000  60.0000  61.0000  62.0000  63.0000  64.0000  65.0000  66.0000
## Effort   163.2622 167.9268 172.5915 177.2561 181.9207 186.5854 191.2500 195.9146 200.5793 205.2439 209.9085 214.5732 219.2378 223.9024 228.5671 233.2317 237.8963 242.5610 247.2256 251.8902 256.5549 261.2195 265.8841 270.5488 275.2134 279.8780 284.5427 289.2073 293.8720 298.5366 303.2012 307.8659
## Richness  81.9161  82.3425  82.7161  83.0964  83.4936  83.8933  84.2293  84.5283  84.8701  85.2724  85.5387  85.8263  86.1494  86.4410  86.7592  87.0225  87.3165  87.6131  87.8547  88.0760  88.3363  88.5027  88.6978  88.9875  89.2456  89.4841  89.7261  89.8954  90.0932  90.3161  90.4846  90.6542
## sd         3.1399   3.0622   3.0889   2.9418   2.8768   2.8643   2.8771   2.9532   2.9370   2.9479   2.9663   2.9171   2.8620   2.8931   2.8170   2.7562   2.7547   2.8049   2.7477   2.7540   2.7294   2.7151   2.7060   2.7082   2.6848   2.6323   2.6701   2.6262   2.6682   2.6914   2.6771   2.6663
##                                                                                                                                                                                                                                                                                                         
## Sites     67.0000  68.0000  69.0000  70.0000  71.0000  72.0000  73.0000  74.0000  75.0000  76.0000  77.0000  78.0000  79.0000  80.0000  81.0000  82.0000  83.0000  84.0000  85.0000  86.0000  87.0000  88.0000  89.0000  90.0000  91.0000  92.0000  93.0000  94.0000  95.0000  96.0000  97.0000  98.0000
## Effort   312.5305 317.1951 321.8598 326.5244 331.1890 335.8537 340.5183 345.1829 349.8476 354.5122 359.1768 363.8415 368.5061 373.1707 377.8354 382.5000 387.1646 391.8293 396.4939 401.1585 405.8232 410.4878 415.1524 419.8171 424.4817 429.1463 433.8110 438.4756 443.1402 447.8049 452.4695 457.1341
## Richness  90.8779  91.0779  91.2097  91.3421  91.5238  91.6833  91.8987  92.1100  92.2622  92.3996  92.4958  92.6370  92.8291  92.9901  93.1603  93.3885  93.5503  93.6618  93.8185  93.9687  94.1282  94.3118  94.4901  94.6310  94.8387  94.9822  95.1005  95.2411  95.3443  95.4566  95.5830  95.7483
## sd         2.6270   2.6534   2.6440   2.5938   2.6117   2.6135   2.6498   2.6679   2.6908   2.6421   2.6006   2.5911   2.6056   2.6737   2.6933   2.5987   2.5340   2.5391   2.5614   2.5676   2.5888   2.5099   2.5361   2.4771   2.4544   2.4693   2.4538   2.4307   2.4065   2.3788   2.3770   2.4033
##                                                                                                                                                                                                                                                                                                        
## Sites     99.0000 100.0000 101.0000 102.0000 103.0000 104.0000 105.0000 106.0000 107.0000 108.0000 109.0000 110.0000 111.0000 112.000 113.0000 114.0000 115.0000 116.0000 117.0000 118.0000 119.0000 120.0000 121.0000 122.0000 123.0000 124.0000 125.0000 126.0000 127.0000 128.0000 129.0000 130.0000
## Effort   461.7988 466.4634 471.1280 475.7927 480.4573 485.1220 489.7866 494.4512 499.1159 503.7805 508.4451 513.1098 517.7744 522.439 527.1037 531.7683 536.4329 541.0976 545.7622 550.4268 555.0915 559.7561 564.4207 569.0854 573.7500 578.4146 583.0793 587.7439 592.4085 597.0732 601.7378 606.4024
## Richness  95.9083  96.0591  96.1962  96.4262  96.6301  96.7518  96.8709  97.0113  97.1618  97.3120  97.4559  97.5820  97.6625  97.774  97.8875  98.0212  98.1916  98.3062  98.4148  98.5660  98.6982  98.8264  98.9860  99.0901  99.2508  99.4181  99.5826  99.6636  99.7552  99.8501  99.9673 100.0626
## sd         2.3710   2.3421   2.3244   2.2608   2.1510   2.1341   2.1339   2.1036   2.0830   2.1229   2.1097   2.0828   2.0669   2.041   2.0307   2.0210   2.0420   2.0676   2.0742   2.0749   2.0260   1.9965   2.0042   1.9885   1.9345   1.9464   1.8630   1.8210   1.7939   1.7916   1.7548   1.7520
##                                                                                                                                                                                                                                                                                                         
## Sites    131.0000 132.0000 133.0000 134.0000 135.0000 136.0000 137.0000 138.0000 139.0000 140.0000 141.0000 142.0000 143.0000 144.0000 145.0000 146.0000 147.0000 148.0000 149.0000 150.0000 151.0000 152.0000 153.0000 154.0000 155.0000 156.0000 157.0000 158.0000 159.0000 160.0000 161.0000 162.0000
## Effort   611.0671 615.7317 620.3963 625.0610 629.7256 634.3902 639.0549 643.7195 648.3841 653.0488 657.7134 662.3780 667.0427 671.7073 676.3720 681.0366 685.7012 690.3659 695.0305 699.6951 704.3598 709.0244 713.6890 718.3537 723.0183 727.6829 732.3476 737.0122 741.6768 746.3415 751.0061 755.6707
## Richness 100.1148 100.2055 100.3376 100.4429 100.5246 100.6333 100.7329 100.7909 100.8665 100.9216 100.9800 101.0768 101.1577 101.2530 101.3564 101.4664 101.5507 101.6486 101.7269 101.7910 101.8817 101.9670 102.0370 102.1180 102.2777 102.3719 102.4443 102.5209 102.5763 102.6302 102.6951 102.7833
## sd         1.7401   1.7160   1.6178   1.6107   1.6278   1.6061   1.6110   1.6030   1.5834   1.5584   1.5354   1.5105   1.4802   1.4457   1.3989   1.3571   1.2741   1.1569   1.1216   1.1036   1.0892   1.0635   1.0139   0.9504   0.8409   0.7648   0.7066   0.6452   0.6152   0.5658   0.5135   0.4589
##                      
## Sites    163.0000 164
## Effort   760.3354 765
## Richness 102.8792 103
## sd         0.3024   0
vegecomm_richness
##          vibifq vibifqnotrees carex cyper dicot shade shrub hydro svp         ap       fqi        bryo   hydro_rel hydro_rel_notrees    sens_rel sens_rel_notree     tol_rel tol_rel_notree         inv  small_tree  subcanopy     canopy  S surveycount
## 1002 47.2642730    24.3506371     9     9    30    25     1     9   1 0.05882353 19.470626 0.005272844 0.015818531       0.031686275 0.039624636     0.000313725 0.419138856    0.766065359 0.000000000 0.139364303 0.03419060 0.19225398 22           3
## 1003 76.1688941    52.6722718     4     4    22    22     1     2   1 0.13333333 22.980970 0.009033832 0.018067663       0.066770532 0.488459280     0.051747162 0.110965566    0.042844425 0.000000000 0.048458150 0.03497367 0.20545653 22           6
## 1004 37.8842226    35.1032304     2     2    21    23     1     9   4 0.19047619 15.086592 0.039605806 0.426344850       0.558710387 0.013978520     0.009159187 0.198727955    0.260426207 0.000000000 0.108695652 0.46386133 0.21743769 28           7
## 1005 82.9198687    65.8474332    11    11    47    52     2    18   3 0.08000000 30.512293 0.007461016 0.060135790       0.128541811 0.168544356     0.088512041 0.133204009    0.196852905 0.000000000 0.114361702 0.10777783 0.17574889 30           3
## 1006 72.0717642    64.7405493     3     3    21    23     1     6   2 0.06666667 26.582281 0.008306874 0.110075316       0.254709325 0.259857491     0.046004015 0.123864727    0.068771091 0.000000000 0.132352941 0.22866471 0.17078115 26           7
## 1008 56.3031565    25.1357383     4     4    11    13     0     3   0 0.25000000 16.842857 0.010954302 0.054771512       0.258754528 0.325945265     0.051750906 0.060577292    0.286182508 0.000000000 0.130841121 0.04519758 0.29052172 24           3
## 1009 45.4587252    35.2832254     5     5    36    40     3    10   0 0.05882353 23.372545 0.005486333 0.109214601       0.114944889 0.010972666     0.003849460 0.611835849    0.643860752 0.000000000 0.119873817 0.12339044 0.13820327 25           3
## 1011 13.3041598    14.0109307     6     6    21    14     0    13   0 0.06976744  9.373602 0.005230308 0.020921232       0.023174523 0.005230308     0.005793631 0.755936399    0.831443966 0.117681928 0.162162162 0.00000000 0.24848942 21           3
## 1012 77.4977507    60.9815207     7     7    28    30     1     6   0 0.04000000 25.285714 0.008523696 0.353818616       0.600318164 0.200562564     0.051050906 0.052505967    0.088796761 0.000000000 0.109965636 0.21892384 0.19566287 21           3
## 1013 48.7604648    40.4112843     7     7    15    22     1     3   0 0.04761905 20.329320 0.031235358 0.078400750       0.100126321 0.078713103     0.040688784 0.258784945    0.329698823 0.000000000 0.219512195 0.00000000 0.27540253 14           3
## 1014 24.8850825    21.4048119     4     4    24    23     1    13   1 0.07142857 12.964988 0.002847867 0.076892408       0.157857811 0.000057000     0.000000000 0.440023922    0.488248363 0.000000000 0.071428571 0.00000000 0.35521137 28           6
## 1015 29.1872412    14.1206164     2     2    13     9     0     2   0 0.10000000 12.159200 0.035009919 0.035009919       0.079702444 0.257206208     0.000000000 0.319990664    0.486184910 0.000000000 0.285714286 0.00000000 0.17139260 26           8
## 1016  6.1773014     6.2905703     2     2    24    12     0    13   2 0.21212121  9.043629 0.000000000 0.015707654       0.016215624 0.006703124     0.006919896 0.902659465    0.931850553 0.001709297 0.071428571 0.00000000 0.27012951 44           3
## 1017 57.8356213    34.9038026     1     2    33    27     1    10   0 0.13793103 22.115051 0.096606419 0.111209715       0.225168877 0.252884151     0.136465986 0.219173004    0.239065663 0.000000000 0.157575758 0.03637554 0.20339806 16           3
## 1018 58.2695721    59.3467604     2     2    38    39     1    12   0 0.10256410 24.999623 0.006024546 0.217838023       0.288496100 0.035789380     0.035548534 0.372955159    0.410980547 0.000000000 0.400000000 0.00000000 0.12057408 33           7
## 1019 50.7262363    18.4208465     4     4    18    16     0     6   0 0.13333333 17.709696 0.009402768 0.018993590       0.074653526 0.018993590     0.036957191 0.181708483    0.455497382 0.000000000 0.290322581 0.03699858 0.15836353 32           7
## 1021 49.5166347    40.4169203     5     5    30    31     1    11   4 0.11764706 18.166431 0.004370948 0.409412108       0.575619149 0.006687550     0.000000000 0.181787718    0.255587193 0.000000000 0.090909091 0.46426353 0.18603507 20           3
## 1022 34.0289308    14.9683193     2     2    17    13     0     1   0 0.08333333 14.058212 0.003540110 0.000000000       0.000000000 0.129618918     0.098373506 0.204447118    0.456594160 0.000000000 0.127819549 0.12032384 0.24034522 35           6
## 1024 54.0397256    25.9199503     2     2    26    20     1     5   0 0.00000000 20.080463 0.003951346 0.007902692       0.022591637 0.248934783     0.011295819 0.216513991    0.404333829 0.000000000 0.081081081 0.00000000 0.26441864 25           3
## 1025 64.6754065    45.8101430     4     4    29    33     1     5   0 0.12500000 23.574758 0.009263548 0.034028099       0.086935942 0.139323761     0.071946986 0.152169214    0.365099400 0.000000000 0.157894737 0.10546062 0.18502495 27           7
## 1026 41.4042709    26.7446869     0     0    20    13     1     6   1 0.22222222 14.971083 0.008514512 0.064903074       0.206543967 0.064260469     0.000000000 0.207079362    0.298057260 0.000000000 0.195555556 0.00000000 0.18063116 27           7
## 1028 50.7125368    32.5246403     5     5    23    24     1    10   3 0.04545454 20.396691 0.006966861 0.094439671       0.314004067 0.092984371     0.007875840 0.479521299    0.428280956 0.000000000 0.161616162 0.09295112 0.19941292 18           3
## 1029  0.8577480     0.8595756     0     0    10     8     0     6   0 0.10000000  6.548462 0.000000000 0.014648816       0.014721654 0.000000000     0.000000000 0.975455957    0.975431413 0.941016686 0.190476190 0.00000000 0.13077558 26           3
## 1030  6.6366366     8.4582943     0     0    12     9     0     5   0 0.33333333  8.249579 0.040040040 0.040840841       0.095595127 0.000000000     0.000000000 0.422422422    0.613870665 0.000000000 0.342105263 0.00000000 0.19200706 25           6
## 1031 86.1822440    76.9382667     5     5    36    41     1     3   3 0.02631579 29.552883 0.009376319 0.016455439       0.030452889 0.291728524     0.515934988 0.070900595    0.113858523 0.000000000 0.139130435 0.03919509 0.14176939 21           3
## 1033 90.5469103    87.9129472     1     1    41    40     2     5   3 0.09375000 31.895261 0.009911785 0.158786798       0.252515434 0.294908646     0.421069224 0.020616513    0.017023512 0.000000000 0.031578947 0.16319081 0.26290521 25           7
## 1034 69.7163721    64.8746123     4     4    55    59     1    23   2 0.07692308 29.634900 0.004247006 0.202030069       0.231641995 0.019238937     0.017189326 0.327444152    0.341351772 0.004247006 0.179104478 0.35823996 0.25887917 28           3
## 1035 69.3713786    45.5529917     9    10    23    30     1    14   2 0.02857143 22.248595 0.006892511 0.034600407       0.085434690 0.427542475     0.042717345 0.131922666    0.257665579 0.000000000 0.230303030 0.11035408 0.19782846 21           3
## 1038 66.9286154    28.3226130     0     0    10    11     1     2   2 0.00000000 16.977494 0.028012512 0.056025025       0.214477212 0.653625286     0.107238606 0.022830197    0.085254692 0.000000000 0.275000000 0.00000000 0.28353301 21           3
## 1039 66.3433378    56.3244555     6     6    55    52     1    27   4 0.08064516 26.131183 0.003231696 0.254323741       0.346897844 0.145426342     0.013224112 0.287987246    0.388406862 0.000000000 0.068181818 0.11838253 0.16168238 22           3
## 1040 34.3195861    25.5902373     2     2    18    14     1    12   1 0.15789474 13.926212 0.011650259 0.648531096       0.713675214 0.020679210     0.000000000 0.255257179    0.280897436 0.000000000 0.324324324 0.00000000 0.22559473 22           3
## 1041 72.4508917    52.4095266     4     4    40    41     1    10   7 0.07500000 26.464671 0.051294008 0.037491257       0.090835483 0.137840989     0.109740712 0.145962540    0.240439113 0.000000000 0.041666667 0.10538593 0.20442155 28           7
## 1043 41.4032009    37.9442588     2     2    18    21     0    14   0 0.13636364 14.234383 0.005054079 0.126402507       0.130355468 0.494036187     0.501667883 0.306024462    0.315594704 0.222379460 0.230769231 0.00000000 0.28416918 30           3
## 1044 19.8247978    11.6533864     0     0     9     5     0     3   1 0.33333333  9.250000 0.026954178 0.026954178       0.099601594 0.000539084     0.000000000 0.567115903    0.601593625 0.000000000 0.276190476 0.00000000 0.20367641 19           3
## 1045 71.7057443    60.6749723    10    10    43    53     1    10   5 0.05660377 27.236578 0.013985199 0.028389954       0.042681431 0.117440709     0.034534114 0.130831537    0.186179346 0.000000000 0.224000000 0.11129807 0.28997505 29           3
## 1046 16.7298468    16.7308176     0     0     1     0     0     5   0 0.20000000  6.123724 0.000000000 0.000000000       0.000000000 0.000000000     0.000000000 0.735939192    0.735981896 0.000000000 0.000000000 0.00000000 0.00000000 32           6
## 1047 91.0710699    71.2182695     5     5    45    48     1    10   3 0.11363636 32.173359 0.015856740 0.025709061       0.078284942 0.381978287     0.126118586 0.038669302    0.098435589 0.000000000 0.134883721 0.12313322 0.21540168 18           3
## 1048 61.6108655    37.9946548     4     4    20    24     1     5   3 0.11764706 19.750829 0.008026326 0.052171121       0.120664583 0.544612997     0.238420841 0.194116703    0.395563256 0.000000000 0.061818182 0.13551587 0.20659002 28           3
## 1050 68.2351091    48.2027445     5     5    25    25     1     8   4 0.10526316 23.859095 0.017453022 0.079062191       0.195877775 0.127494328     0.087345056 0.095060795    0.192274431 0.000000000 0.150943396 0.03500834 0.17121292 21           7
## 1051 60.0000000    25.0000000     0     0     5     6     1     2   0 0.00000000 14.000000 0.010141580 0.030123506       0.298804781 0.692840647     0.000000000 0.000000000    0.000000000 0.000000000 0.107954545 0.00000000 0.21661381 27           7
## 1052 21.8469209    19.8484848     0     0     7     5     0     1   0 0.50000000  6.717514 0.000000000 0.005478852       0.009090909 0.000000000     0.000000000 0.139710717    0.231818182 0.000000000 0.135338346 0.00000000 0.40077160 23           6
## 1053 79.1463360    75.8082659     2     2    54    54     1    20   1 0.08620690 31.544663 0.006684045 0.171735401       0.295427542 0.215861239     0.083649534 0.073234855    0.125982139 0.000000000 0.280000000 0.11111865 0.11976083 33           7
## 1054 49.3529079    42.0231102     0     0    25    22     1     8   0 0.10000000 18.007351 0.009740971 0.280661737       0.345487430 0.017046700     0.002997722 0.023183512    0.025540589 0.000000000 0.197674419 0.07054649 0.28838175 32           7
## 1055  5.7153765     4.5936048     1     1     6     6     0     2   0 0.12500000  4.472136 0.002974951 0.000000000       0.000000000 0.000000000     0.000000000 0.845005057    0.977156566 0.000000000 0.046321526 0.00000000 0.24992981 27           6
## 1057 47.0121070    39.0228034     1     1    18    25     1     5   4 0.10526316 18.762424 0.091991619 0.209843103       0.251269812 0.061327746     0.055076189 0.296826289    0.355425005 0.000000000 0.380952381 0.00000000 0.15870464 22           3
## 1058 32.2451969    26.7484609     1     1    20    16     1     9   0 0.20000000 13.139781 0.006768266 0.151157937       0.412142835 0.020304798     0.000000000 0.091913050    0.250238366 0.000000000 0.133333333 0.27578122 0.26959703 15           3
## 1059  7.4596278     7.7087769     4     4    16     8     0     6   0 0.16666667  5.208554 0.000000000 0.004988327       0.005544221 0.000000000     0.000000000 0.819921384    0.911292469 0.000000000 0.000000000 0.00000000 0.12568772 26           3
## 1060 86.0217979    75.2021735     2     2    27    33     1     3   5 0.00000000 29.340949 0.008283327 0.193277620       0.317710067 0.128962191     0.189294684 0.011375769    0.018699507 0.000000000 0.074380165 0.31247163 0.26409609 22           3
## 1061 64.4979646    44.4512500     7     8    24    29     1     8   1 0.03571429 23.590713 0.021376149 0.078379213       0.160279761 0.214331520     0.088008160 0.185806237    0.234251299 0.000000000 0.145038168 0.05658729 0.18345411 19           7
## 1063  0.6882786     0.7016067     0     0     7     3     0     1   0 0.23076923  2.910427 0.003754247 0.000000000       0.000000000 0.000000000     0.000000000 0.969665684    0.988442625 0.000000000 0.000000000 0.00000000 0.00000000 38           7
## 1064 68.9815761    48.7614881     3     3    14    15     1     2   2 0.09090909 20.199962 0.008518371 0.014197285       0.129735340 0.259810324     0.233523612 0.000851837    0.007784120 0.000000000 0.082901554 0.03687149 0.20108903 34           8
## 1065 36.4339155    26.4460066     6     7    41    25     2    30   3 0.04000000 19.398969 0.005679667 0.166641424       0.194026232 0.005679667     0.002204343 0.745708696    0.808736544 0.000000000 0.056547619 0.00000000 0.17096747 25           3
## 1066 71.6540458    68.8886354     4     4    29    31     1     2   3 0.00000000 27.955103 0.034103698 0.038650858       0.098776909 0.097229643     0.178757154 0.051360169    0.104935940 0.000000000 0.140625000 0.11569459 0.25493198 25           3
## 1067 49.0975392    28.8159047     2     2    19    18     1     5   0 0.10000000 17.826752 0.011034686 0.116084894       0.402602373 0.176996359     0.000000000 0.227424872    0.213930348 0.000000000 0.054347826 0.14472548 0.16772444 22           7
## 1068 86.8148642    77.5296486     6     6    51    53     1    17   5 0.14583333 30.137684 0.049954209 0.134412503       0.163033599 0.385931943     0.394362133 0.096958740    0.108948743 0.000000000 0.241379310 0.08966525 0.20454119 23           8
## 1069 49.5071980    51.5128576     5     5    41    32     1    21   1 0.11904762 22.584706 0.002664960 0.200584444       0.310688341 0.072692731     0.106464500 0.204456553    0.298294729 0.000000000 0.152777778 0.23370100 0.12741797 40           3
## 1070 69.8973526    55.2878470     2     2    20    20     2    10   5 0.00000000 23.334524 0.068522032 0.311195446       0.557770430 0.248875536     0.085025980 0.132932743    0.030420406 0.000000000 0.132911392 0.30080543 0.17851444 14           3
## 1071 44.0556446    31.0731281     3     3    28    24     1    11   1 0.09523809 18.577086 0.017223067 0.076585240       0.157366993 0.202371042     0.026542409 0.185664667    0.381148991 0.000000000 0.148148148 0.03828634 0.09913688 30           7
## 1072 22.7373661    21.3626378     5     6    17    17     1     8   0 0.18421053 12.450604 0.011351382 0.281230490       0.383464465 0.000681083     0.000619115 0.394668634    0.390790662 0.000000000 0.340425532 0.00000000 0.39070309 31           7
## 1073 32.0698544     7.3519025     6     6    35    16     1     7   0 0.09375000 15.748536 0.002367527 0.004782405       0.008101442 0.101803661     0.000000000 0.534212738    0.884908725 0.000000000 0.030470914 0.00000000 0.19053064 35           7
## 1074 30.4423741    13.7202290     0     0    26    15     0     1   0 0.20000000 14.411534 0.004714313 0.000000000       0.000000000 0.009428625     0.000000000 0.410522346    0.673015873 0.000000000 0.162500000 0.00000000 0.42606512 46           7
## 1075 30.9225792    19.1232646     4     5    26    20     0    12   1 0.03225807 14.142857 0.024070607 0.064348756       0.151197134 0.024070607     0.056557532 0.546563252    0.624018098 0.000000000 0.242424242 0.00000000 0.30032618 33           6
## 1076 35.1904406    29.1257065     5     6    34    22     1    22   2 0.07142857 17.008401 0.009766367 0.377361562       0.406555154 0.015878485     0.003507305 0.630220937    0.652671544 0.000000000 0.103734440 0.30960972 0.20621328 25           3
## 1077 49.8635649    24.4490187     8     8    27    21     1    11   2 0.12000000 17.649664 0.004873769 0.159640316       0.441313183 0.063359002     0.013473155 0.289136368    0.704982822 0.000000000 0.092807425 0.17718386 0.20787458 17           3
## 1078 26.2286053    22.1832371     4     4    26    24     1    15   0 0.21875000 14.010778 0.002774310 0.107273316       0.116382900 0.124843945     0.135445616 0.561011698    0.608652467 0.123456790 0.210526316 0.08238879 0.08245940 29           3
## 1080 22.1493130    14.1714373     5     6    28    20     0    14   1 0.03333333 15.100260 0.003637488 0.032737392       0.058499805 0.006529291     0.000000000 0.418602113    0.631017897 0.000000000 0.125954198 0.03460696 0.16864491 32           3
## 1081 67.6441016    33.1461639     4     4    23    21     1     8   0 0.33333333 19.776610 0.006473051 0.120830277       0.520156047 0.514866439     0.028422813 0.181633798    0.558424670 0.000000000 0.067524116 0.06885338 0.47382135 26           7
## 1082 52.2306642    42.5636909     5     5    43    40     1    13   4 0.04444444 24.002841 0.018320051 0.061097371       0.107320998 0.022961131     0.016090105 0.440464922    0.467882363 0.000000000 0.099630996 0.00000000 0.27425889 33           7
## 1083 32.2417471    25.9868553     1     2    13    13     0     9   0 0.07142857 12.701706 0.008325124 0.308029582       0.665228335 0.004245813     0.000000000 0.142359618    0.217188062 0.000000000 0.161290323 0.00000000 0.19785790 29           7
## 1086 59.1996968    52.8073944     1     1    32    26     1    18   0 0.03333333 21.144229 0.004861291 0.636018927       0.708483755 0.048612912     0.043321300 0.112619912    0.114620939 0.009722582 0.196721311 0.47819559 0.17963810 34           6
## 1087 69.9907734    38.0826786     2     2    13    13     0     2   1 0.10000000 18.766739 0.010051261 0.000201025       0.001755926 0.563071665     0.089552239 0.007370925    0.064383963 0.000000000 0.132530120 0.00000000 0.32294425 18           3
## 1088 22.6745597    23.3575929     0     0    11     3     0     2   0 0.07692308 10.253048 0.003005651 0.000000000       0.000000000 0.026048972     0.026856175 0.439947101    0.453580135 0.000000000 0.000000000 0.00000000 0.00000000 25           3
## 1089 92.2211128    61.0635030     4     4    44    41     1     9   6 0.19444444 30.250000 0.087474376 0.029423199       0.088148874 0.603467166     0.183286126 0.069493532    0.096699050 0.000000000 0.095238095 0.00000000 0.28637885 26           7
## 1090 55.8924083    31.2538117     3     3    34    29     1     9   0 0.17391304 20.164033 0.017208150 0.027808370       0.095081196 0.323788546     0.047540598 0.175201909    0.551973013 0.000000000 0.208695652 0.22057155 0.17443767 37           7
## 1091 55.3715779    33.3228587     6     6    13    15     1     6   1 0.00000000 14.818657 0.015372134 0.208411973       0.718504343 0.389529869     0.017665244 0.031359153    0.107757986 0.000000000 0.130434783 0.14917110 0.20109531 19           3
## 1092 63.1114138    19.1674678     3     3    15    12     0     1   0 0.00000000 18.557783 0.019445457 0.007058242       0.042725913 0.342395332     0.064516129 0.077922995    0.471694082 0.000000000 0.156950673 0.03491555 0.18949034 19           3
## 1093 56.1844921    26.1134860     3     3    19    13     0     4   1 0.09090909 18.278154 0.009517766 0.019225888       0.179396092 0.095748731     0.092362345 0.191878173    0.454706927 0.000000000 0.272727273 0.00000000 0.15406053 23           3
## 1094 34.6268174    21.0755361     1     1    33    24     1     9   1 0.10000000 17.259516 0.010561521 0.031790178       0.081041059 0.018641084     0.000000000 0.393416652    0.437514023 0.000000000 0.086419753 0.03497019 0.17006012 31           3
## 1097 75.3499787    45.2044752     2     2    27    21     0     4   1 0.10526316 24.833333 0.005980980 0.017942941       0.131262306 0.495015850     0.055348939 0.041986483    0.307153796 0.000000000 0.078804348 0.13952624 0.16486691 33           7
## 1098 65.0847028    62.3960785     4     4    58    50     1    21   2 0.06349206 26.052335 0.002918827 0.171354628       0.190698933 0.003093957     0.003378269 0.240355708    0.228444589 0.002918827 0.162790698 0.03799072 0.13694672 41           7
## 1099 57.6102830    20.0327567     3     3    15    11     0     1   0 0.12500000 16.341190 0.010840892 0.000000000       0.000000000 0.553319120     0.057164207 0.156614751    0.417336073 0.000000000 0.201550388 0.03598872 0.22731842 24           3
## 1100 28.3592783     7.1167580     1     1    12     7     0     5   1 0.00000000 11.783766 0.014668492 0.029336984       0.074738416 0.132309799     0.000747384 0.514668492    0.861235675 0.000000000 0.236714976 0.00000000 0.19891336 20           3
## 1101 66.6230627    53.7345943     5     5    43    39     0     0   2 0.04761905 27.287642 0.002503703 0.000000000       0.000000000 0.047057105     0.025676038 0.309232406    0.358469679 0.000000000 0.013392857 0.17047783 0.20196864 24           7
## 1102 43.5536105    26.8481837     5     5    21    27     0     0   0 0.05555556 16.302936 0.011794926 0.000000000       0.000000000 0.009304886     0.024861685 0.179544978    0.479305273 0.007863284 0.220512821 0.00000000 0.14265512 33           6
## 1103 51.0794746    23.7739794     8     9    22    23     0     0   2 0.06250000 18.243355 0.007918910 0.000000000       0.000000000 0.079347482     0.042271382 0.219934537    0.587008595 0.000000000 0.131034483 0.03510817 0.21502618 12           1
## 1104 17.4910984    14.0366974     6     6    31    22     0    20   1 0.08571429 15.105187 0.002750313 0.027503128       0.029980063 0.005500626     0.002998006 0.841664489    0.863500772 0.000000000 0.060498221 0.00000000 0.27189190 30           3
## 1105 64.6964124    54.1651226     8     8    27    33     0     0   1 0.07692308 24.651184 0.015704340 0.000000000       0.000000000 0.171988693     0.170685279 0.206040936    0.177203507 0.000000000 0.222222222 0.06094992 0.18350302 29           7
## 1107  7.5175044     6.5512624     0     0    12     9     0     0   0 0.94736842  4.585303 0.000000000 0.000000000       0.000000000 0.018754689     0.000000000 0.772693173    0.788064269 0.018754689 0.000000000 0.00000000 0.00000000 25           3
## 1108 39.3594051    36.4378180     2     2    34    22     2    39   1 0.13953488 20.579116 0.000000000 0.101119398       0.103930836 0.072136131     0.070169864 0.775862735    0.797434166 0.000000000 0.000000000 0.16789173 0.30721036 32           3
## 1110 50.3361915    38.9448141     3     3    16    17     0     0   0 0.07692308 16.614251 0.012599483 0.000000000       0.000000000 0.021083136     0.000000000 0.056277693    0.138987289 0.000000000 0.116666667 0.42389991 0.21260174 20           3
## 1111 40.3099379    36.7915692     2     2    21    23     0     0   1 0.05555556 18.584678 0.019406171 0.000000000       0.000000000 0.077624685     0.079978672 0.349311081    0.399893362 0.000000000 0.233766234 0.00000000 0.21955821 19           3
## 1112 49.5745936    26.5228283     6     6    30    28     1    15   2 0.05405405 19.564659 0.032777482 0.098496333       0.198590153 0.180276150     0.049564930 0.286283990    0.461559643 0.000000000 0.123076923 0.19252234 0.21024158 26           3
## 1113 67.3125523    44.0752952     1     1    19    17     1     1   0 0.11111111 22.299904 0.006195275 0.004171485       0.038055765 0.260325458     0.169555388 0.241739633    0.170685757 0.000000000 0.052380952 0.04775089 0.24970445 27           7
## 1114 38.5917371    14.6986197     2     2    16    14     0     0   1 0.16666667 15.493175 0.012769354 0.000000000       0.000000000 0.221744736     0.000000000 0.241696851    0.469809069 0.000000000 0.203883495 0.00000000 0.17589600 21           7
## 1115 37.0430682    16.6196300     1     1    11    10     0     0   0 0.10000000 13.764944 0.033090668 0.000000000       0.000000000 0.083057578     0.047977422 0.036399735    0.101599247 0.000000000 0.068493151 0.00000000 0.22685374 17           3
## 1116 22.0782247    22.8754990     5     5    27    21     0     0   2 0.18181818 13.567731 0.003262536 0.000000000       0.000000000 0.009787609     0.010380623 0.613356824    0.650519031 0.000000000 0.000000000 0.00000000 0.00000000 33           6
## 1117 80.3722481    54.2961888     3     3    24    30     0     0   2 0.14814815 24.363107 0.008905116 0.000000000       0.000000000 0.445612004     0.105351171 0.023420455    0.109113712 0.000000000 0.047732697 0.03492895 0.21164164 24           7
## 1119  5.0552703     4.7242690     2     2    15     4     1     3   0 0.07692308  8.966120 0.000000000 0.000000000       0.000000000 0.011494693     0.000000000 0.889516840    0.968385919 0.000000000 0.018867925 0.00000000 0.15788103 27           6
## 1122 49.1094343    36.7583160     9    10    22    25     0     0   0 0.02631579 20.178224 0.012632430 0.000000000       0.000000000 0.071667986     0.051813770 0.322564889    0.441015780 0.012632430 0.383561644 0.00000000 0.18456142 30           6
## 1123  3.9409783     3.9054676     4     5    16    10     0     0   1 0.12500000  8.198916 0.000000000 0.000000000       0.000000000 0.001956954     0.001963698 0.974999918    0.978360051 0.013698675 0.000000000 0.00000000 0.00000000 36           3
## 1124 26.7138853     6.3973468     3     4    26    19     0    10   2 0.00000000 14.320780 0.005453703 0.016579256       0.031289133 0.223710881     0.000205850 0.485470427    0.878394373 0.000000000 0.079787234 0.03600028 0.11806617 28           3
## 1125 44.0553690    32.2007370     7     7    30    32     0     0   3 0.09375000 20.445485 0.012575242 0.000000000       0.000000000 0.041917473     0.000000000 0.448282222    0.557006820 0.000000000 0.136842105 0.03690905 0.14177715 26           3
## 1126 26.0212404    16.7881897     4     4    23    17     0     0   1 0.10000000 14.142857 0.017190992 0.000000000       0.000000000 0.086642599     0.045945325 0.530256146    0.593728463 0.000000000 0.315217391 0.00000000 0.13719234 32           6
## 1127 68.5741717    41.9298087     1     1     8    11     0     0   2 0.00000000 18.432707 0.016758279 0.000000000       0.000000000 0.643718997     0.139013008 0.015384100    0.045576408 0.000000000 0.118644068 0.00000000 0.43130675 37           8
## 1128 64.7408687    52.9120284     2     2    23    25     1     3   2 0.09090909 23.509055 0.019483682 0.055528495       0.085675635 0.272771554     0.210431384 0.166975158    0.257628138 0.000000000 0.063492063 0.13589695 0.26069283 19           6
## 1129 74.6561865    66.0798092     0     0    15    18     1     1   1 0.09090909 23.337820 0.011290717 0.000225814       0.000634350 0.367625750     0.667336258 0.045388683    0.064069356 0.000000000 0.056074766 0.08701678 0.20152191 27           7
## 1131 45.5803365    25.6452461     1     1    12    10     0     0   3 0.09090909 15.652476 0.026589852 0.000000000       0.000000000 0.020607135     0.000000000 0.080833149    0.188819876 0.000000000 0.081395349 0.00000000 0.29849496 28           3
## 1132 61.5461263    48.7898138     0     0    29    32     0     0   3 0.03571429 23.193102 0.008254912 0.000000000       0.000000000 0.181608057     0.068073519 0.227587915    0.301452235 0.000000000 0.205882353 0.44181765 0.28576944 23           6
## 1133 11.0775751    11.0812355     2     3    17     5     0     0   0 0.07407407  6.761234 0.000000000 0.000000000       0.000000000 0.001251951     0.001255961 0.778054134    0.780520889 0.000000000 0.000000000 0.00000000 0.00000000 31           3
## 1134  6.3143417     4.1146735     3     3    17    12     0     0   0 0.10526316  9.328719 0.010534633 0.000000000       0.000000000 0.010534633     0.000000000 0.758581336    0.867277354 0.000000000 0.224719101 0.00000000 0.20337584 31           7
## 1135 62.4428095    56.3290470     7     7    37    34     0     0   0 0.05263158 22.648834 0.005368983 0.000000000       0.000000000 0.025270015     0.012427073 0.162375954    0.149145382 0.000000000 0.108391608 0.24237264 0.18391901 30           7
## 1136 48.4135472    26.0122411     1     1    10     4     0     0   0 0.00000000 15.588457 0.018674718 0.000000000       0.000000000 0.467241433     0.141242938 0.299168975    0.002824859 0.000000000 0.312169312 0.03520657 0.17819956 36           7
## 1138 14.4219204    11.3714328     3     3    35    20     0     0   0 0.07317073 11.324417 0.002655854 0.000000000       0.000000000 0.002655854     0.002834226 0.704624728    0.751948530 0.000000000 0.003267974 0.00000000 0.20428325 38           7
## 1140 48.2773792    35.8349200     3     3    11    15     1     1   2 0.07692308 17.800000 0.024457851 0.024457851       0.041006014 0.305723137     0.164024057 0.082830589    0.138873701 0.000000000 0.207207207 0.00000000 0.30972068  9           2
## 1141 19.2394574    19.7753194     7     9    12    15     0     0   2 0.07142857 11.833333 0.005878780 0.000000000       0.000000000 0.017636339     0.018281164 0.832337207    0.862769393 0.005878780 0.416243655 0.07767635 0.06882081 35           4
## 1143 56.2833929    44.5067290     3     3    21    26     0     0   2 0.04761905 22.062113 0.013690191 0.000000000       0.000000000 0.082414950     0.059209881 0.113697036    0.162827174 0.000000000 0.157894737 0.03578353 0.22902443 24           3
## 1144 64.4036676    50.7355473     7     7    29    36     1    12   3 0.02631579 24.270696 0.009465215 0.132513015       0.243013366 0.104306673     0.052074293 0.223284430    0.270612741 0.000000000 0.046920821 0.12178217 0.15508403 18           2
## 1145 74.0871563    55.5889607     2     2    29    29     1     9   3 0.07692308 24.675454 0.027309516 0.143052211       0.376281591 0.227438616     0.082348331 0.062563619    0.164566055 0.000000000 0.013245033 0.08775583 0.30325680 27           7
## 1146 48.7429834    33.8597926     7     7    23    25     0     0   2 0.00000000 20.833333 0.042099354 0.000000000       0.000000000 0.042099354     0.000000000 0.151557676    0.179257362 0.000000000 0.285714286 0.00000000 0.17181709 20           3
## 1147 55.8073915    35.4732695     3     3    12    14     0     0   1 0.11111111 18.371173 0.022140221 0.000000000       0.000000000 0.067306273     0.057475582 0.045166052    0.114951165 0.000000000 0.094339623 0.07925455 0.15954906 26           7
## 1148 31.3509544    11.9440054     0     0    19    10     0     3   0 0.00000000 12.247449 0.012371644 0.000000000       0.000000000 0.024743288     0.000000000 0.435399398    0.644840822 0.000000000 0.162162162 0.14737364 0.23031636 20           3
## 1149 31.6438090    14.7058823     2     2    15    12     0     0   0 0.14285714 13.761951 0.018796992 0.000000000       0.000000000 0.125563910     0.000000000 0.206766917    0.294117647 0.000000000 0.171171171 0.00000000 0.23571960 39           7
## 1150 13.7657362     9.2470760     0     0     7     5     0     0   0 0.00000000  6.948792 0.082101806 0.000000000       0.000000000 0.082101806     0.000000000 0.415435140    0.554824561 0.000000000 0.562500000 0.00000000 0.17073555 29           6
## 1151 53.0986317    19.0950517     3     3    12    10     0     0   1 0.00000000 15.600270 0.012645955 0.000000000       0.000000000 0.089027526     0.004713276 0.038443704    0.358208955 0.000000000 0.118421053 0.00000000 0.17254395 20           3
## 1153 64.9326107    51.9177502     2     2    27    29     0     0   3 0.12500000 24.019223 0.020253848 0.000000000       0.000000000 0.175533351     0.149812734 0.108290575    0.200249688 0.000000000 0.103448276 0.00000000 0.42040109 29           7
## 1154 55.5187565    45.7347820     0     0    21    21     0     0   0 0.05263158 16.537377 0.004964504 0.000000000       0.000000000 0.084396565     0.021809062 0.047659236    0.052341748 0.000000000 0.022421525 0.36813310 0.13428338 28           3
## 1155 40.3122142    24.3452125     2     2    12    15     0     0   1 0.07142857 12.894156 0.012861185 0.000000000       0.000000000 0.285175341     0.000000000 0.142244705    0.237440962 0.000000000 0.137254902 0.00000000 0.23696839 24           3
## 1156  1.0813491     1.0265884     1     1    17    14     0    12   0 0.09523809  9.128709 0.002620923 0.015725537       0.015892145 0.000000000     0.000000000 0.952770971    0.962865354 0.002620923 0.478260870 0.00000000 0.24308641 30           3
## 1157 61.5372140    40.1781105     7     7    15    13     0     0   0 0.11764706 19.843135 0.010522439 0.000000000       0.000000000 0.226653338     0.142510121 0.081583978    0.192442645 0.000000000 0.071770335 0.03514090 0.19232775 17           3
## 1159 34.1045005    20.7453541     6     6    23    22     0     0   1 0.06896552 15.776277 0.029364551 0.000000000       0.000000000 0.047218198     0.010911252 0.367879096    0.466861125 0.000000000 0.217777778 0.10655868 0.19307106 16           1
## 1161 44.2238010    36.0447104     2     2    14    14     1     7   3 0.10000000 16.819419 0.027068483 0.135883786       0.262827225 0.198502211     0.104712042 0.231976902    0.291623037 0.000000000 0.240000000 0.00000000 0.19752813 29           7
## 1162 84.7094202    82.8188118     7     8    61    49     0     0   4 0.09090909 28.781504 0.006421576 0.000000000       0.000000000 0.404719859     0.477399404 0.284582865    0.281107633 0.000000000 0.166666667 0.11502740 0.16604261 27           3
## 1164 27.9469239    10.6513497     2     2    23    16     0     7   0 0.00000000 12.809880 0.009602612 0.016260423       0.025736505 0.177936399     0.000000000 0.418401805    0.601134838 0.009602612 0.061538462 0.00000000 0.20613724 24           3
## 1166 46.4687101    42.2716973     2     2    20    22     0     0   2 0.04166667 19.438548 0.021276596 0.000000000       0.000000000 0.063829787     0.051948052 0.212765957    0.259740260 0.021276596 0.269230769 0.00000000 0.26124182 28           3
## 1167 15.6939345     9.4909865     4     5    24    17     0     0   1 0.02564103 13.483997 0.001710888 0.000000000       0.000000000 0.005166883     0.003714768 0.861229844    0.934935843 0.000000000 0.055016181 0.08762413 0.17993393 20           3
## 1168 33.4930442    32.1376096     2     2    25    25     1    12   0 0.02941176 16.130857 0.002667876 0.262172185       0.265014428 0.005335752     0.005393598 0.539617960    0.545468029 0.002667876 0.046511628 0.08721168 0.21476748 25           3
## 1169 44.8892920    35.6223573     5     5    24    27     0     0   2 0.11111111 19.975541 0.016398371 0.000000000       0.000000000 0.049523081     0.045172219 0.311897018    0.339243365 0.000000000 0.277777778 0.00000000 0.23198093 20           3
## 1171 51.1157162    32.4117032     2     2    21    26     0     0   1 0.13636364 18.493242 0.007276167 0.000000000       0.000000000 0.109288027     0.094694479 0.146360098    0.317463240 0.000000000 0.164473684 0.18532007 0.21561005 24           3
## 1173 11.7999514     4.9824524     1     1    32    21     0     0   3 0.20408163 13.076697 0.002379272 0.000000000       0.000000000 0.009517087     0.000000000 0.873704288    0.927745079 0.130859948 0.120879121 0.00000000 0.17192061 27           3
## 1174 45.4563394    41.3060154     6     8    38    31     0     0   2 0.10204082 17.094791 0.004708172 0.000000000       0.000000000 0.009416344     0.005153311 0.355466972    0.353001804 0.000000000 0.063492063 0.55001752 0.21682028 22           3
## 1175 47.1780885    30.1299671     2     2    26    24     0     0   0 0.20689655 19.990242 0.020109934 0.000000000       0.000000000 0.127631050     0.087834870 0.228448854    0.374176548 0.000000000 0.082568807 0.03527536 0.20087728 18           3
## 1176 79.3295904    62.5540727     5     5    22    26     1     3   4 0.05555556 27.666667 0.023475097 0.017723698       0.040644206 0.361105677     0.134852631 0.039281662    0.049706159 0.000000000 0.082758621 0.03568503 0.18198608 25           6
## 1177 10.8313493     4.2607287     1     1    21    16     0     7   1 0.00000000 11.672328 0.002994610 0.013475744       0.015087170 0.020962268     0.000000000 0.659612697    0.735136343 0.000000000 0.036659878 0.00000000 0.24035969 21           3
## 1178 63.5212313    48.0175070     5     5    15    20     0     0   3 0.00000000 23.186922 0.017584994 0.000000000       0.000000000 0.211019930     0.100446429 0.070691676    0.067633929 0.000000000 0.094936709 0.12220624 0.18210034 29           7
## 1179 19.2360846    11.6142806     3     3    24    19     0     0   0 0.10810811 13.152186 0.000000000 0.000000000       0.000000000 0.009867123     0.008879242 0.470727536    0.552959747 0.002466781 0.093525180 0.00000000 0.25779149 26           3
## 1180  1.8083853     2.5876619     4     4    13    13     0     0   0 0.00000000  9.449112 0.003971564 0.000000000       0.000000000 0.003971564     0.005689253 0.644386195    0.923081299 0.003971564 0.062015504 0.00000000 0.25432460 30           6
## 1182 47.4456174    32.5749792     0     0    16    13     0     0   0 0.11111111 17.000000 0.019193858 0.000000000       0.000000000 0.038771593     0.063511831 0.192706334    0.250311333 0.000000000 0.283185841 0.00000000 0.18594614 32           7
## 1183 22.1725624    15.4762436     2     2    16    11     1     4   0 0.00000000 10.253048 0.014928344 0.044785032       0.080303368 0.014928344     0.000000000 0.293665406    0.526031675 0.000000000 0.212765957 0.00000000 0.26602949 23           6
## 1184 36.0002961    20.7327971     6     6    10     7     0     6   0 0.16666667 12.578842 0.028764562 0.014669927       0.024859859 0.273263340     0.000000000 0.179778513    0.255910310 0.000000000 0.105485232 0.00000000 0.23273561 16           3
## 1187 31.4765682    24.5520339     2     2    28    24     0     0   2 0.14285714 17.195706 0.008974915 0.000000000       0.000000000 0.049362033     0.042464562 0.511570161    0.655167533 0.000000000 0.017953321 0.00000000 0.20523612 22           3
## 1188 11.4124737    11.6851910     1     1    22    11     0     9   1 0.06451613  9.170605 0.003078328 0.328416559       0.343207639 0.001067154     0.001115216 0.918732139    0.934373844 0.000000000 0.020408163 0.00000000 0.09499185 27           3
## 1189 78.2998025    69.5002046    12    12    49    59     0     0   7 0.08571429 31.045780 0.006506632 0.000000000       0.000000000 0.176980402     0.089774282 0.233436283    0.265128035 0.000000000 0.188888889 0.15676817 0.13104936 30           3
## 1190 23.0916312    21.6842615     5     5    26    29     0     0   0 0.10256410 16.424201 0.004528097 0.000000000       0.000000000 0.031696678     0.032531408 0.754682807    0.774557340 0.566012105 0.166666667 0.05848294 0.12946399 32           6
## 1191 48.9535176    37.1858220     5     5    24    26     0     0   2 0.04166667 20.739806 0.012896570 0.000000000       0.000000000 0.077379417     0.046468401 0.202390164    0.364622057 0.000000000 0.059090909 0.07114911 0.11167694 28           8
## 1192 72.4031045    50.2791237     4     5    37    32     1    12   5 0.07894737 23.576676 0.014490654 0.294643288       0.532611543 0.403323190     0.034925347 0.103994590    0.187985681 0.000000000 0.014150943 0.22526683 0.32548705 30           6
## 1194 21.3756974     8.5813381     2     2    27    17     0     0   2 0.05555556 14.388363 0.006841505 0.000000000       0.000000000 0.006841505     0.000000000 0.464082098    0.564371257 0.000000000 0.065868263 0.08676337 0.22470780 21           3
## 1195 24.4356468     6.7525974     3     4    30    24     0     0   0 0.10526316 14.293825 0.004343231 0.000000000       0.000000000 0.078178159     0.000000000 0.616912542    0.829813786 0.000000000 0.047272727 0.00000000 0.20777354 27           7
## 1196 33.4535675    23.5771189     6     6    34    29     0     0   0 0.15384615 16.016063 0.006009435 0.000000000       0.000000000 0.186292479     0.054861657 0.497631281    0.656017409 0.000000000 0.144067797 0.00000000 0.22926804 35           8
## 1197 74.0828500    58.7189028     1     1    27    34     0     0   4 0.07142857 26.231569 0.020425764 0.000000000       0.000000000 0.272479688     0.092109303 0.082247742    0.139085048 0.000000000 0.049180328 0.25528860 0.26281393 24           3
## 1198 10.5315333     9.5710090     3     4    19    11     0     0   0 0.00000000 10.224154 0.012729124 0.000000000       0.000000000 0.000000000     0.000000000 0.681262729    0.763373191 0.000000000 0.129921260 0.00000000 0.27265831 20           3
## 1199 72.7472391    68.0281154     8     8    27    32     0     0   2 0.09677419 24.917691 0.042283960 0.000000000       0.000000000 0.260218624     0.279594541 0.156654242    0.179143236 0.000000000 0.323809524 0.03651024 0.13588334 19           3
## 1200 51.7404400    39.5929323     1     1     9     9     0     0   0 0.25000000 15.377079 0.130737943 0.000000000       0.000000000 0.581057525     0.489457831 0.022225450    0.001129518 0.000000000 0.155555556 0.18603820 0.23530338 32           6
## 1257 62.8582229    36.3095238     1     1    14    12     1     1   2 0.09090909 21.059956 0.025000000 0.025000000       0.047619048 0.437500000     0.190476190 0.125000000    0.238095238 0.000000000 0.188888889 0.00000000 0.36842214 32           3
## 3980  1.1376107     1.0261715     0     0     4     3     0     0   0 0.00000000  2.474874 0.000000000 0.000000000       0.000000000 0.000000000     0.000000000 0.951658235    0.960496503 0.000000000 0.086956522 0.00000000 0.21715732 23           6
CMP_bbs_pcap_join_pcaponly_pcapsum_date<-CMP_bbs_pcap_join_pcaponly %>% 
  group_by(simp_plot, reservatio,date,species_4_) %>%
  summarise(n=n(),longitude=mean(xcoord), latitude=mean(ycoord),vibifq=mean(vibi_fq),vibifqnotrees=mean(vibi_fq_no_trees),carex=mean(carex_metric_value),cyper=mean(cyperaceae_metric_value),
            dicot=mean(dicot_metric_value),shade=mean(shade_metric_value),shrub=mean(shrub_metric_value),
            hydro=mean(hydrophyte_metric_value),svp=mean(svp_metric_value),ap=mean(ap_ratio_metric_value),
            fqi=mean(fqai_score),bryo=mean(bryophyte_metric_value),hydro_rel=mean(hydrophyte_rel_cov_metric_value),
            hydro_rel_notrees=mean(hydrophyte_rel_cov_no_trees),
            sens_rel=mean(sensitive_rel_cov_metric_value),
            sens_rel_notree=mean(sensitive_rel_cov_no_trees),
            tol_rel=mean(tolerant_rel_cov_metric_value),
            tol_rel_notree=mean(tolerant_rel_cov_no_trees),
            inv=mean(invasive_graminoids_metric_value),
            small_tree=mean(small_tree_metric_value),
            subcanopy=mean(subcanopy_iv),
            canopy=mean(canopy_iv.1)
            ) %>%
  spread(species_4_,n) %>%
  replace(is.na(.),0)

birds_bycount<-CMP_bbs_pcap_join_pcaponly_pcapsum_date[,29:ncol(CMP_bbs_pcap_join_pcaponly_pcapsum_date)]

vegecomm_bycount<-CMP_bbs_pcap_join_pcaponly_pcapsum_date[,1:28]

index<-as.numeric(factor(CMP_bbs_pcap_join_pcaponly_pcapsum_date$simp_plot))
index_unique<-unique(index)

rarefylist<-list()

for(i in index_unique){
  birds_temp<-birds_bycount[index %in% i,]
  rarefylist[[i]]<- specaccum(birds_temp)
}

totalspecaccum<-specaccum(birds_bycount)

plot(totalspecaccum)

plot(rarefylist[[1]])

So the number of surveys clearly has an effect on the species richness of the site. This is somewhat difficult to account for, but we will attempt to do so using mixed effects models. I then use Akaike’s Information Criterion (AIC for short) to rank these models. The lowest number is the best, and any model within 2 AIC of the top model is potentially worth interpreting.

simp_plot<-factor(CMP_bbs_pcap_join_pcaponly_pcapsum_date$simp_plot)

date<-strptime(CMP_bbs_pcap_join_pcaponly_pcapsum_date$date,format="%m/%d/%Y")

month<-month(date)

S_plotdate<-specnumber(birds_bycount)


Rich_vibifq<-lmer(S_plotdate~vibifq+(1|simp_plot),vegecomm_bycount)
summary(Rich_vibifq)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ vibifq + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4350.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6215 -0.6926 -0.1170  0.6499  3.4453 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  5.692   2.386   
##  Residual              13.710   3.703   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 14.42728    0.51881  27.809
## vibifq      -0.02997    0.01009  -2.972
## 
## Correlation of Fixed Effects:
##        (Intr)
## vibifq -0.892
Rich_vibifqnotrees<-lmer(S_plotdate~vibifqnotrees+(1|simp_plot),vegecomm_bycount)
summary(Rich_vibifqnotrees)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ vibifqnotrees + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4356.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5584 -0.6948 -0.1026  0.6376  3.4659 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.088   2.467   
##  Residual              13.702   3.702   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)   13.61824    0.46208   29.47
## vibifqnotrees -0.01692    0.01192   -1.42
## 
## Correlation of Fixed Effects:
##             (Intr)
## vibifqnotrs -0.855
Rich_carex<-lmer(S_plotdate~carex+(1|simp_plot),vegecomm_bycount)
summary(Rich_carex)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ carex + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4353.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5078 -0.6891 -0.0956  0.6328  3.5149 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.196   2.489   
##  Residual              13.685   3.699   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 12.74937    0.40130  31.770
## carex        0.09391    0.09746   0.964
## 
## Correlation of Fixed Effects:
##       (Intr)
## carex -0.800
Rich_cyper<-lmer(S_plotdate~cyper+(1|simp_plot),vegecomm_bycount)
summary(Rich_cyper)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ cyper + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4353.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5164 -0.6870 -0.0952  0.6341  3.5192 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.186   2.487   
##  Residual              13.685   3.699   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 12.71752    0.39948   31.84
## cyper        0.09949    0.09302    1.07
## 
## Correlation of Fixed Effects:
##       (Intr)
## cyper -0.798
Rich_dicot<-lmer(S_plotdate~dicot+(1|simp_plot),vegecomm_bycount)
summary(Rich_dicot)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ dicot + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4351.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4704 -0.7081 -0.0930  0.6339  3.5457 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  5.934   2.436   
##  Residual              13.667   3.697   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 11.80814    0.54551  21.646
## dicot        0.05250    0.02068   2.539
## 
## Correlation of Fixed Effects:
##       (Intr)
## dicot -0.900
Rich_shade<-lmer(S_plotdate~shade+(1|simp_plot),vegecomm_bycount)
summary(Rich_shade)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ shade + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4356.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4711 -0.6893 -0.0932  0.6487  3.5347 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.176   2.485   
##  Residual              13.676   3.698   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 12.48609    0.49794  25.075
## shade        0.02613    0.01990   1.313
## 
## Correlation of Fixed Effects:
##       (Intr)
## shade -0.876
Rich_hydro<-lmer(S_plotdate~hydro+(1|simp_plot),vegecomm_bycount)
summary(Rich_hydro)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ hydro + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4349.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4429 -0.6824 -0.0808  0.6376  3.4771 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  5.831   2.415   
##  Residual              13.684   3.699   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 12.52655    0.30454  41.132
## hydro        0.09721    0.03549   2.739
## 
## Correlation of Fixed Effects:
##       (Intr)
## hydro -0.632
Rich_shrub<-lmer(S_plotdate~shrub+(1|simp_plot),vegecomm_bycount)
summary(Rich_shrub)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ shrub + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4351.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5380 -0.6891 -0.0877  0.6446  3.4829 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.204   2.491   
##  Residual              13.691   3.700   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  13.2033     0.3119   42.34
## shrub        -0.3018     0.4136   -0.73
## 
## Correlation of Fixed Effects:
##       (Intr)
## shrub -0.635
Rich_svp<-lmer(S_plotdate~svp+(1|simp_plot),vegecomm_bycount)
summary(Rich_svp)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ svp + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4353.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5198 -0.6806 -0.1005  0.6550  3.4904 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.248   2.500   
##  Residual              13.686   3.699   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 13.02335    0.31742  41.029
## svp          0.02676    0.15316   0.175
## 
## Correlation of Fixed Effects:
##     (Intr)
## svp -0.649
Rich_ap<-lmer(S_plotdate~ap+(1|simp_plot),vegecomm_bycount)
summary(Rich_ap)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ ap + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4348.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5238 -0.6858 -0.0978  0.6528  3.4990 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.234   2.497   
##  Residual              13.691   3.700   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  13.0950     0.3456  37.896
## ap           -0.3580     2.4652  -0.145
## 
## Correlation of Fixed Effects:
##    (Intr)
## ap -0.716
Rich_bryo<-lmer(S_plotdate~bryo+(1|simp_plot),vegecomm_bycount)
summary(Rich_bryo)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ bryo + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4333.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5734 -0.6798 -0.0879  0.6431  3.5273 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  5.611   2.369   
##  Residual              13.679   3.699   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  13.6941     0.2986  45.862
## bryo        -42.5099    12.3301  -3.448
## 
## Correlation of Fixed Effects:
##      (Intr)
## bryo -0.625
Rich_sens_rel<-lmer(S_plotdate~sens_rel+(1|simp_plot),vegecomm_bycount)
summary(Rich_sens_rel)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ sens_rel + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4337.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5731 -0.7099 -0.1163  0.6427  3.5445 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  5.461   2.337   
##  Residual              13.724   3.705   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  13.7862     0.3107  44.375
## sens_rel     -4.8222     1.3559  -3.556
## 
## Correlation of Fixed Effects:
##          (Intr)
## sens_rel -0.668
Rich_tol_rel<-lmer(S_plotdate~tol_rel+(1|simp_plot),vegecomm_bycount)
summary(Rich_tol_rel)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ tol_rel + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4330.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6961 -0.6952 -0.1092  0.6514  3.4089 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  5.137   2.267   
##  Residual              13.698   3.701   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  11.8052     0.3544  33.310
## tol_rel       3.8867     0.8554   4.544
## 
## Correlation of Fixed Effects:
##         (Intr)
## tol_rel -0.769
Rich_fqi<-lmer(S_plotdate~fqi+(1|simp_plot),vegecomm_bycount)
summary(Rich_fqi)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ fqi + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4355.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5737 -0.6954 -0.1016  0.6440  3.4662 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.134   2.477   
##  Residual              13.698   3.701   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 13.83929    0.70282  19.691
## fqi         -0.04354    0.03680  -1.183
## 
## Correlation of Fixed Effects:
##     (Intr)
## fqi -0.940
Rich_inv<-lmer(S_plotdate~inv+(1|simp_plot),vegecomm_bycount)
summary(Rich_inv)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ inv + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4346.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5113 -0.6818 -0.0954  0.6528  3.5191 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.115   2.473   
##  Residual              13.697   3.701   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  13.0095     0.2425  53.658
## inv           3.7372     2.8028   1.333
## 
## Correlation of Fixed Effects:
##     (Intr)
## inv -0.148
Rich_small_tree<-lmer(S_plotdate~small_tree+(1|simp_plot),vegecomm_bycount)
summary(Rich_small_tree)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ small_tree + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4348.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5315 -0.6775 -0.0946  0.6508  3.4892 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.242   2.498   
##  Residual              13.688   3.700   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  13.1300     0.4135  31.750
## small_tree   -0.4937     2.3432  -0.211
## 
## Correlation of Fixed Effects:
##            (Intr)
## small_tree -0.812
Rich_subcanopy<-lmer(S_plotdate~subcanopy+(1|simp_plot),vegecomm_bycount)
summary(Rich_subcanopy)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ subcanopy + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4348.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5203 -0.6845 -0.0962  0.6450  3.4966 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.238   2.498   
##  Residual              13.689   3.700   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  13.0349     0.2907   44.85
## subcanopy     0.3244     2.1587    0.15
## 
## Correlation of Fixed Effects:
##           (Intr)
## subcanopy -0.557
Rich_canopy<-lmer(S_plotdate~canopy+(1|simp_plot),vegecomm_bycount)
summary(Rich_canopy)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ canopy + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4340.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7009 -0.6938 -0.0975  0.6317  3.6263 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  5.823   2.413   
##  Residual              13.686   3.699   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  14.7055     0.6461  22.762
## canopy       -7.9969     2.9116  -2.747
## 
## Correlation of Fixed Effects:
##        (Intr)
## canopy -0.931
Rich_null<-lmer(S_plotdate~1+(1|simp_plot),vegecomm_bycount)
summary(Rich_null)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ 1 + (1 | simp_plot)
##    Data: vegecomm_bycount
## 
## REML criterion at convergence: 4351.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5202 -0.6845 -0.0981  0.6522  3.5041 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  simp_plot (Intercept)  6.176   2.485   
##  Residual              13.691   3.700   
## Number of obs: 765, groups:  simp_plot, 164
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  13.0584     0.2406   54.28
richmods<-list(Rich_ap,Rich_bryo,Rich_canopy,Rich_carex,Rich_cyper,Rich_dicot,Rich_fqi,Rich_hydro,Rich_inv,Rich_sens_rel,Rich_shade,Rich_shrub,Rich_small_tree,Rich_subcanopy,Rich_svp,Rich_tol_rel,Rich_vibifq,Rich_vibifqnotrees,Rich_null)

aictab(richmods)
## 
## Model selection based on AICc:
## 
##       K    AICc Delta_AICc AICcWt Cum.Wt   Res.LL
## Mod16 4 4338.87       0.00   0.76   0.76 -2165.41
## Mod2  4 4341.52       2.65   0.20   0.96 -2166.73
## Mod10 4 4345.35       6.48   0.03   0.99 -2168.65
## Mod3  4 4348.56       9.69   0.01   1.00 -2170.26
## Mod9  4 4354.29      15.43   0.00   1.00 -2173.12
## Mod1  4 4356.31      17.44   0.00   1.00 -2174.13
## Mod13 4 4356.39      17.52   0.00   1.00 -2174.17
## Mod14 4 4356.57      17.70   0.00   1.00 -2174.26
## Mod8  4 4357.41      18.54   0.00   1.00 -2174.68
## Mod19 3 4357.95      19.08   0.00   1.00 -2175.96
## Mod17 4 4358.70      19.83   0.00   1.00 -2175.33
## Mod12 4 4359.36      20.50   0.00   1.00 -2175.66
## Mod6  4 4359.50      20.63   0.00   1.00 -2175.72
## Mod5  4 4361.74      22.87   0.00   1.00 -2176.84
## Mod15 4 4361.86      22.99   0.00   1.00 -2176.90
## Mod4  4 4361.86      22.99   0.00   1.00 -2176.90
## Mod7  4 4363.34      24.47   0.00   1.00 -2177.64
## Mod11 4 4364.24      25.37   0.00   1.00 -2178.09
## Mod18 4 4364.98      26.11   0.00   1.00 -2178.46

So richness - at a plot level, and after accounting for plot identity - is determined by how many tolerant species are in the associated plot. And this is much more predictive than any other model.

Here is the plot for this top model:

## `geom_smooth()` using formula 'y ~ x'

We can also verify that our mixed effects model fits better compared to the fixed effects model.

Rich_tol_rel_fixed<-lm(S_plotdate~tol_rel,vegecomm_bycount)
summary(Rich_tol_rel_fixed)
## 
## Call:
## lm(formula = S_plotdate ~ tol_rel, data = vegecomm_bycount)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.515  -2.830  -0.563   2.757  19.963 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.6653     0.2401  48.585  < 2e-16 ***
## tol_rel       3.7349     0.5971   6.255  6.6e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.29 on 763 degrees of freedom
## Multiple R-squared:  0.04878,    Adjusted R-squared:  0.04754 
## F-statistic: 39.13 on 1 and 763 DF,  p-value: 6.604e-10
anova(Rich_tol_rel, Rich_tol_rel_fixed) 
## refitting model(s) with ML (instead of REML)
## Data: vegecomm_bycount
## Models:
## Rich_tol_rel_fixed: S_plotdate ~ tol_rel
## Rich_tol_rel: S_plotdate ~ tol_rel + (1 | simp_plot)
##                    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)    
## Rich_tol_rel_fixed    3 4403.2 4417.1 -2198.6   4397.2                        
## Rich_tol_rel          4 4339.2 4357.8 -2165.6   4331.2 66.02  1  4.464e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This confirms that we chose correctly - the mixed effects model fits better than the fixed effects-only model.

We can also inspect how this observation-level, mixed effect model compares to a more consolidated, plot-level view (in other words - a view that just looks at richness at a plot level - aggregating over observations).

S_tol_rel<-lm(S~tol_rel,vegecomm_richness)
confint(S_tol_rel)
##                 2.5 %   97.5 %
## (Intercept) 22.889198 25.86972
## tol_rel      1.996856  9.10019
confint(Rich_tol_rel)
## Computing profile confidence intervals ...
##                 2.5 %    97.5 %
## .sig01       1.832966  2.688447
## .sigma       3.500278  3.923018
## (Intercept) 11.111432 12.501114
## tol_rel      2.211719  5.563805

We can see that, while the covariate estimate for the plot-level model is higher, it is also estimated with much more variance (its confidence intervals are wider). This is not surprising - our functional sample size for this dataset is 164 plots. In contrast, while the observation-level mixed-effects model has a lower covariate estimate, we have tighter confidence intervals around this estimate - likely by virtue of a larger number of samples and a more effectively parameterized model.

reservation-level analyses

birdcomm_reserv<-CMP_bbs_pcap_join_pcaponly %>%
  group_by(reservatio,species_4_) %>%
  summarise(n=n()) %>%
  spread(species_4_,n) %>%
  replace(is.na(.),0) 
## `summarise()` has grouped output by 'reservatio'. You can override using the `.groups` argument.
birdcomm_reserv<-column_to_rownames(birdcomm_reserv,var="reservatio")
birdcomm_reserv<-birdcomm_reserv[,-1]
birdcomm_reserv_bin<-birdcomm_reserv %>% mutate_if(is.numeric, ~1 * (. >= 1))

reservspcount<-specnumber(birdcomm_reserv_bin)
reservspcount
## BC BE BK BR BW EC GP HI HU MS NC OE RR SC WA WC 
## 48 73 45 70 41 20 24 78 27 63 66 43 80 76 23 48
vegecomm_reserv<-CMP_bbs_pcap_join_pcaponly %>%
  group_by(reservatio) %>%
  summarise(vibifq=mean(vibi_fq),vibifqnotrees=mean(vibi_fq_no_trees),carex=mean(carex_metric_value),cyper=mean(cyperaceae_metric_value),
            dicot=mean(dicot_metric_value),shade=mean(shade_metric_value),shrub=mean(shrub_metric_value),
            hydro=mean(hydrophyte_metric_value),svp=mean(svp_metric_value),ap=mean(ap_ratio_metric_value),
            fqi=mean(fqai_score),bryo=mean(bryophyte_metric_value),hydro_rel=mean(hydrophyte_rel_cov_metric_value),
            hydro_rel_notrees=mean(hydrophyte_rel_cov_no_trees),
            sens_rel=mean(sensitive_rel_cov_metric_value),
            sens_rel_notree=mean(sensitive_rel_cov_no_trees),
            tol_rel=mean(tolerant_rel_cov_metric_value),
            tol_rel_notree=mean(tolerant_rel_cov_no_trees),
            inv=mean(invasive_graminoids_metric_value),
            small_tree=mean(small_tree_metric_value),
            subcanopy=mean(subcanopy_iv),
            canopy=mean(canopy_iv.1),
            numsurveys=n_distinct(paste0(simp_plot,date)))
vegecomm_reserv_richness<-cbind(vegecomm_reserv,reservspcount)

cor(vegecomm_reserv_richness[,-1])
##                        vibifq vibifqnotrees       carex       cyper        dicot       shade      shrub       hydro         svp           ap          fqi        bryo   hydro_rel hydro_rel_notrees    sens_rel sens_rel_notree      tol_rel tol_rel_notree         inv  small_tree    subcanopy
## vibifq             1.00000000     0.9241914  0.50266385  0.46630888  0.564528005  0.65220135  0.7527538  0.62273545  0.60880098  0.257794968  0.947059034  0.30558082  0.62339581        0.65183940  0.56704683      0.42599572 -0.852896083     -0.8105305  0.32149298  0.47496212  0.531367247
## vibifqnotrees      0.92419139     1.0000000  0.53929085  0.47464154  0.449543250  0.61107030  0.7782038  0.55620231  0.69989878  0.418705980  0.874506913  0.36484997  0.59042482        0.60034564  0.52764476      0.53330180 -0.774869734     -0.8047408  0.28706233  0.40312084  0.498283056
## carex              0.50266385     0.5392908  1.00000000  0.98861968  0.572023490  0.75776669  0.5042649  0.35192082  0.64080449  0.377693834  0.671921171 -0.25412469  0.32776255        0.34946703  0.07325813     -0.08076780 -0.304555329     -0.2892609  0.53268817  0.18937253 -0.281897425
## cyper              0.46630888     0.4746415  0.98861968  1.00000000  0.584049479  0.75372315  0.4711865  0.30703125  0.58053080  0.380730464  0.647650588 -0.27287406  0.31028593        0.32971714  0.06376714     -0.10310014 -0.258492814     -0.2380678  0.50561201  0.17547331 -0.327850311
## dicot              0.56452800     0.4495432  0.57202349  0.58404948  1.000000000  0.94021548  0.4698399  0.54034866  0.48543907  0.008029017  0.759772038 -0.24608019  0.28671560        0.24788748 -0.08333476     -0.19172516 -0.347413659     -0.1282112  0.17164936 -0.03846365  0.155163341
## shade              0.65220135     0.6110703  0.75776669  0.75372315  0.940215480  1.00000000  0.5829774  0.52014187  0.66723418  0.142508957  0.832756899 -0.22197887  0.35664035        0.33764329 -0.01115955     -0.10193730 -0.369664790     -0.2287267  0.35303648  0.06045936  0.084086166
## shrub              0.75275378     0.7782038  0.50426492  0.47118653  0.469839900  0.58297741  1.0000000  0.78513894  0.61806333  0.111257132  0.750745575 -0.13100838  0.83342517        0.83262041  0.23501494      0.14386524 -0.460370277     -0.4995103  0.10638055  0.41413586  0.249525499
## hydro              0.62273545     0.5562023  0.35192082  0.30703125  0.540348664  0.52014187  0.7851389  1.00000000  0.40397499 -0.209004300  0.622184752 -0.17324938  0.62466824        0.59390990  0.20199468      0.01283236 -0.409560618     -0.3488116  0.12178203  0.20959270  0.183422757
## svp                0.60880098     0.6998988  0.64080449  0.58053080  0.485439069  0.66723418  0.6180633  0.40397499  1.00000000  0.124468114  0.688362107 -0.04239540  0.33225582        0.38825400  0.11004921      0.13158407 -0.319130764     -0.3035986  0.33680986  0.06588272  0.216164439
## ap                 0.25779497     0.4187060  0.37769383  0.38073046  0.008029017  0.14250896  0.1112571 -0.20900430  0.12446811  1.000000000  0.268493012  0.46101347  0.16767876        0.14567773  0.42055027      0.56317362 -0.425818603     -0.4889413  0.04450250  0.09192957  0.037297463
## fqi                0.94705903     0.8745069  0.67192117  0.64765059  0.759772038  0.83275690  0.7507456  0.62218475  0.68836211  0.268493012  1.000000000  0.12111170  0.57120898        0.58946433  0.39402289      0.24845247 -0.761357994     -0.6636032  0.30819955  0.33480585  0.391252116
## bryo               0.30558082     0.3648500 -0.25412469 -0.27287406 -0.246080192 -0.22197887 -0.1310084 -0.17324938 -0.04239540  0.461013466  0.121111704  1.00000000 -0.11350270       -0.10247675  0.76072190      0.92221351 -0.529298731     -0.5976709 -0.02147223  0.02772700  0.472984957
## hydro_rel          0.62339581     0.5904248  0.32776255  0.31028593  0.286715602  0.35664035  0.8334252  0.62466824  0.33225582  0.167678756  0.571208983 -0.11350270  1.00000000        0.97600979  0.34652164      0.09194741 -0.511414589     -0.5937322 -0.10008676  0.61633226  0.193693179
## hydro_rel_notrees  0.65183940     0.6003456  0.34946703  0.32971714  0.247887485  0.33764329  0.8326204  0.59390990  0.38825400  0.145677729  0.589464327 -0.10247675  0.97600979        1.00000000  0.39250298      0.10495976 -0.534687160     -0.6144103 -0.01962444  0.58390153  0.178383496
## sens_rel           0.56704683     0.5276448  0.07325813  0.06376714 -0.083334765 -0.01115955  0.2350149  0.20199468  0.11004921  0.420550273  0.394022888  0.76072190  0.34652164        0.39250298  1.00000000      0.82206470 -0.684201135     -0.7837929  0.06059223  0.21381404  0.222296394
## sens_rel_notree    0.42599572     0.5333018 -0.08076780 -0.10310014 -0.191725160 -0.10193730  0.1438652  0.01283236  0.13158407  0.563173621  0.248452466  0.92221351  0.09194741        0.10495976  0.82206470      1.00000000 -0.518781157     -0.6364926  0.01042068  0.05071249  0.383357352
## tol_rel           -0.85289608    -0.7748697 -0.30455533 -0.25849281 -0.347413659 -0.36966479 -0.4603703 -0.40956062 -0.31913076 -0.425818603 -0.761357994 -0.52929873 -0.51141459       -0.53468716 -0.68420113     -0.51878116  1.000000000      0.9366517 -0.18307775 -0.45979912 -0.580636269
## tol_rel_notree    -0.81053054    -0.8047408 -0.28926090 -0.23806784 -0.128211218 -0.22872670 -0.4995103 -0.34881162 -0.30359859 -0.488941346 -0.663603226 -0.59767089 -0.59373217       -0.61441028 -0.78379286     -0.63649263  0.936651732      1.0000000 -0.16840035 -0.58104234 -0.502322014
## inv                0.32149298     0.2870623  0.53268817  0.50561201  0.171649356  0.35303648  0.1063806  0.12178203  0.33680986  0.044502499  0.308199552 -0.02147223 -0.10008676       -0.01962444  0.06059223      0.01042068 -0.183077752     -0.1684003  1.00000000  0.23345191 -0.168146221
## small_tree         0.47496212     0.4031208  0.18937253  0.17547331 -0.038463654  0.06045936  0.4141359  0.20959270  0.06588272  0.091929570  0.334805846  0.02772700  0.61633226        0.58390153  0.21381404      0.05071249 -0.459799124     -0.5810423  0.23345191  1.00000000  0.306801643
## subcanopy          0.53136725     0.4982831 -0.28189742 -0.32785031  0.155163341  0.08408617  0.2495255  0.18342276  0.21616444  0.037297463  0.391252116  0.47298496  0.19369318        0.17838350  0.22229639      0.38335735 -0.580636269     -0.5023220 -0.16814622  0.30680164  1.000000000
## canopy             0.02686002    -0.1477253 -0.16296479 -0.12875828  0.004174685 -0.07217575  0.1982734  0.12798031  0.01062822 -0.321878438 -0.009884374 -0.27391776  0.52889987        0.54700481  0.01841165     -0.24518828  0.006449976     -0.0197229 -0.24910337  0.35704003 -0.006191979
## numsurveys         0.67749638     0.7157262  0.43821784  0.38388893  0.511769451  0.64351866  0.7158910  0.51650992  0.66735153  0.213173431  0.671919649 -0.06541184  0.60459521        0.63984087  0.15097169      0.10308858 -0.480145398     -0.4446043  0.30324470  0.21309096  0.354918723
## reservspcount      0.68233393     0.7884269  0.50420389  0.45376839  0.464086302  0.60821512  0.7523661  0.54918263  0.62534630  0.449063177  0.673416476  0.08143411  0.66313940        0.63780487  0.24936734      0.27304031 -0.548039099     -0.5601424  0.23875891  0.30577390  0.344432018
##                         canopy  numsurveys reservspcount
## vibifq             0.026860021  0.67749638    0.68233393
## vibifqnotrees     -0.147725299  0.71572622    0.78842694
## carex             -0.162964792  0.43821784    0.50420389
## cyper             -0.128758278  0.38388893    0.45376839
## dicot              0.004174685  0.51176945    0.46408630
## shade             -0.072175750  0.64351866    0.60821512
## shrub              0.198273381  0.71589097    0.75236610
## hydro              0.127980314  0.51650992    0.54918263
## svp                0.010628216  0.66735153    0.62534630
## ap                -0.321878438  0.21317343    0.44906318
## fqi               -0.009884374  0.67191965    0.67341648
## bryo              -0.273917763 -0.06541184    0.08143411
## hydro_rel          0.528899874  0.60459521    0.66313940
## hydro_rel_notrees  0.547004812  0.63984087    0.63780487
## sens_rel           0.018411651  0.15097169    0.24936734
## sens_rel_notree   -0.245188281  0.10308858    0.27304031
## tol_rel            0.006449976 -0.48014540   -0.54803910
## tol_rel_notree    -0.019722898 -0.44460427   -0.56014242
## inv               -0.249103369  0.30324470    0.23875891
## small_tree         0.357040034  0.21309096    0.30577390
## subcanopy         -0.006191979  0.35491872    0.34443202
## canopy             1.000000000  0.07617249   -0.01191566
## numsurveys         0.076172489  1.00000000    0.90376963
## reservspcount     -0.011915663  0.90376963    1.00000000
ggplot(aes(x=vibifq,y=reservspcount),data=vegecomm_reserv_richness)+geom_point()+geom_smooth(method="lm")
## `geom_smooth()` using formula 'y ~ x'

reservation_spaccum<-specaccum(birdcomm_reserv_bin,w = vegecomm_reserv$numsurveys, method="random")
plot(reservation_spaccum)

This doesn’t quite get at what we want - this displays species accumulation across reservations after accounting for survey effort per reservation. Let’s try exploring reservation-by-reservation rarefaction with the previous observation-level dataset.

birds_bycount
## # A tibble: 765 x 104
##     ACFL  ALFL  AMCR  AMGO  AMKE  AMRE  AMRO  AMWO  BADO  BAEA  BANS  BAOR  BARS  BCCH  BEKI  BGGN  BHCO  BHVI  BLJA  BOBO  BRCR  BRTH  BTNW  BWHA  BWWA  CANG  CARW  CEDW  CERW  CHSP  CHSW  CLSW  COGR  COHA  COYE  DEJU  DOWO  EABL  EAKI  EAME  EAPH  EATO  EAWP  EUST  FISP  GBHE  GCFL  GHOW  GRCA
##    <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
##  1     0     0     1     1     0     1     0     0     0     0     0     0     0     0     0     1     0     0     1     0     0     0     0     0     1     0     0     1     0     0     0     0     0     0     0     0     1     0     0     0     0     1     1     0     0     0     0     0     0
##  2     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     1     1     0     0     0     0     0     0
##  3     1     0     0     0     0     1     0     0     0     0     0     0     0     1     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0
##  4     0     0     0     1     0     0     1     0     0     0     0     1     0     1     0     0     1     0     1     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
##  5     1     0     0     1     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
##  6     1     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     1     0     0
##  7     1     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0
##  8     1     0     0     0     0     0     0     0     0     0     0     1     0     1     0     1     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     1     0     0     0     0     0     0
##  9     1     0     0     1     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0
## 10     1     0     1     0     0     0     1     0     0     0     0     0     0     1     1     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     1     0     0     0     0     0     1     0     0     0     1     0     0
## # ... with 755 more rows, and 55 more variables: GRHE <int>, HAWO <int>, HERG <int>, HETH <int>, HOFI <int>, HOSP <int>, HOWA <int>, HOWR <int>, INBU <int>, KILL <int>, LOWA <int>, MALL <int>, MODO <int>, NOCA <int>, NOFL <int>, NOPA <int>, NRWS <int>, OROR <int>, OVEN <int>, PIWO <int>,
## #   PUFI <int>, PUMA <int>, RBGR <int>, RBGU <int>, RBNU <int>, RBWO <int>, REVI <int>, RHWO <int>, ROPI <int>, RSHA <int>, RTHA <int>, RTHU <int>, RWBL <int>, SAVS <int>, SCTA <int>, SOSP <int>, SWSP <int>, TEWA <int>, TRES <int>, TUTI <int>, TUVU <int>, VEER <int>, VIRA <int>, WAVI <int>,
## #   WBNU <int>, WEVI <int>, WIFL <int>, WITU <int>, WIWR <int>, WODU <int>, WOTH <int>, YBCU <int>, YEWA <int>, YTVI <int>, YTWA <int>
vegecomm_bycount
## # A tibble: 765 x 28
## # Groups:   simp_plot, reservatio, date [765]
##    simp_plot reservatio date      longitude latitude vibifq vibifqnotrees carex cyper dicot shade shrub hydro   svp     ap   fqi    bryo hydro_rel hydro_rel_notrees sens_rel sens_rel_notree tol_rel tol_rel_notree   inv small_tree subcanopy canopy    V1
##        <int> <chr>      <chr>         <dbl>    <dbl>  <dbl>         <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>   <dbl>     <dbl>             <dbl>    <dbl>           <dbl>   <dbl>          <dbl> <dbl>      <dbl>     <dbl>  <dbl> <int>
##  1      1002 HI         6/18/2019  2181968.  561931.   47.3          24.4     9     9    30    25     1     9     1 0.0588  19.5 0.00527    0.0158            0.0317   0.0396        0.000314   0.419         0.766      0     0.139     0.0342  0.192     0
##  2      1002 HI         6/25/2019  2181968.  561931.   47.3          24.4     9     9    30    25     1     9     1 0.0588  19.5 0.00527    0.0158            0.0317   0.0396        0.000314   0.419         0.766      0     0.139     0.0342  0.192     0
##  3      1002 HI         6/5/2019   2181968.  561931.   47.3          24.4     9     9    30    25     1     9     1 0.0588  19.5 0.00527    0.0158            0.0317   0.0396        0.000314   0.419         0.766      0     0.139     0.0342  0.192     0
##  4      1003 MS         6/13/2017  2150599.  613095.   76.2          52.7     4     4    22    22     1     2     1 0.133   23.0 0.00903    0.0181            0.0668   0.488         0.0517     0.111         0.0428     0     0.0485    0.0350  0.205     0
##  5      1003 MS         6/24/2018  2150599.  613095.   76.2          52.7     4     4    22    22     1     2     1 0.133   23.0 0.00903    0.0181            0.0668   0.488         0.0517     0.111         0.0428     0     0.0485    0.0350  0.205     0
##  6      1003 MS         6/25/2017  2150599.  613095.   76.2          52.7     4     4    22    22     1     2     1 0.133   23.0 0.00903    0.0181            0.0668   0.488         0.0517     0.111         0.0428     0     0.0485    0.0350  0.205     0
##  7      1003 MS         6/4/2017   2150599.  613095.   76.2          52.7     4     4    22    22     1     2     1 0.133   23.0 0.00903    0.0181            0.0668   0.488         0.0517     0.111         0.0428     0     0.0485    0.0350  0.205     1
##  8      1003 MS         6/7/2018   2150599.  613095.   76.2          52.7     4     4    22    22     1     2     1 0.133   23.0 0.00903    0.0181            0.0668   0.488         0.0517     0.111         0.0428     0     0.0485    0.0350  0.205     0
##  9      1003 MS         7/1/2018   2150599.  613095.   76.2          52.7     4     4    22    22     1     2     1 0.133   23.0 0.00903    0.0181            0.0668   0.488         0.0517     0.111         0.0428     0     0.0485    0.0350  0.205     0
## 10      1004 BE         6/12/2018  2227458.  626184.   37.9          35.1     2     2    21    23     1     9     4 0.190   15.1 0.0396     0.426             0.559    0.0140        0.00916    0.199         0.260      0     0.109     0.464   0.217     0
## # ... with 755 more rows
index<-as.numeric(factor(CMP_bbs_pcap_join_pcaponly_pcapsum_date$reservatio))
index_unique<-unique(index)

rarefylist_reservation<-list()

for(i in index_unique){
  birds_temp<-birds_bycount[index %in% i,]
  rarefylist_reservation[[i]]<- specaccum(birds_temp)
}
## Warning in cor(x > 0): the standard deviation is zero

## Warning in cor(x > 0): the standard deviation is zero

## Warning in cor(x > 0): the standard deviation is zero

## Warning in cor(x > 0): the standard deviation is zero

## Warning in cor(x > 0): the standard deviation is zero
plot(totalspecaccum)

plot(rarefylist_reservation[[4]])

interpolatedrichness_n3<-data.frame(levels(factor(CMP_bbs_pcap_join_pcaponly_pcapsum_date$reservatio)),rep(1,16))

colnames(interpolatedrichness_n3)<-c("reservation","richness")

for(i in index_unique){
  birds_temp<-birds_bycount[index %in% i,]
  interpolatedrichness_n3$richness[i]<- predict(specaccum(birds_temp),newdata=3)
}
## Warning in cor(x > 0): the standard deviation is zero

## Warning in cor(x > 0): the standard deviation is zero

## Warning in cor(x > 0): the standard deviation is zero

## Warning in cor(x > 0): the standard deviation is zero

## Warning in cor(x > 0): the standard deviation is zero

With the plot-level data, we could easily account for variation in survey effort across plots with a mixed model design, but this case is a bit more complicated. We are interested in reservation-level species richness, but know that reservations were visited different numbers of times - based on the number of plots and the number of surveys/plot. Thus, the easiest thing to start with is simply using species richness estimates interpolated to the least number of surveys per reservation (3, in our case).

Now we can proceed with exploring richness patterns across reservations.

vegecomm_reserv_richness$interrich_n3<-interpolatedrichness_n3$richness

Rich_reserv_vibifq<-lm(interrich_n3~vibifq,vegecomm_reserv_richness)
summary(Rich_reserv_vibifq)
## 
## Call:
## lm(formula = interrich_n3 ~ vibifq, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0917 -1.9946  0.4847  1.9275  4.8225 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 19.24757    1.87384  10.272 6.69e-08 ***
## vibifq       0.12874    0.04718   2.729   0.0163 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.731 on 14 degrees of freedom
## Multiple R-squared:  0.3472, Adjusted R-squared:  0.3005 
## F-statistic: 7.445 on 1 and 14 DF,  p-value: 0.01632
Rich_reserv_vibifqnotrees<-lm(interrich_n3~vibifqnotrees,vegecomm_reserv_richness)
summary(Rich_reserv_vibifqnotrees)
## 
## Call:
## lm(formula = interrich_n3 ~ vibifqnotrees, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3776 -1.4771  0.0605  1.4821  4.7685 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   19.81963    1.40534  14.103 1.15e-09 ***
## vibifqnotrees  0.16551    0.04961   3.336   0.0049 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.523 on 14 degrees of freedom
## Multiple R-squared:  0.4429, Adjusted R-squared:  0.4031 
## F-statistic: 11.13 on 1 and 14 DF,  p-value: 0.004899
Rich_reserv_carex<-lm(interrich_n3~carex,vegecomm_reserv_richness)
summary(Rich_reserv_carex)
## 
## Call:
## lm(formula = interrich_n3 ~ carex, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7041 -2.3742  0.1555  1.4905  5.9512 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.2721     1.6801  13.257 2.58e-09 ***
## carex         0.6488     0.5506   1.178    0.258    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.224 on 14 degrees of freedom
## Multiple R-squared:  0.09023,    Adjusted R-squared:  0.02525 
## F-statistic: 1.388 on 1 and 14 DF,  p-value: 0.2583
Rich_reserv_cyper<-lm(interrich_n3~cyper,vegecomm_reserv_richness)
summary(Rich_reserv_cyper)
## 
## Call:
## lm(formula = interrich_n3 ~ cyper, data = vegecomm_reserv_richness)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.091 -2.482  0.218  1.558  5.838 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.5436     1.7022  13.243 2.61e-09 ***
## cyper         0.5155     0.5252   0.982    0.343    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.269 on 14 degrees of freedom
## Multiple R-squared:  0.06438,    Adjusted R-squared:  -0.00245 
## F-statistic: 0.9633 on 1 and 14 DF,  p-value: 0.343
Rich_reserv_dicot<-lm(interrich_n3~dicot,vegecomm_reserv_richness)
summary(Rich_reserv_dicot)
## 
## Call:
## lm(formula = interrich_n3 ~ dicot, data = vegecomm_reserv_richness)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.896 -1.993  0.058  1.958  4.709 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.9424     2.2412   9.344 2.15e-07 ***
## dicot         0.1489     0.1019   1.461    0.166    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.148 on 14 degrees of freedom
## Multiple R-squared:  0.1324, Adjusted R-squared:  0.07039 
## F-statistic: 2.136 on 1 and 14 DF,  p-value: 0.166
Rich_reserv_shade<-lm(interrich_n3~shade,vegecomm_reserv_richness)
summary(Rich_reserv_shade)
## 
## Call:
## lm(formula = interrich_n3 ~ shade, data = vegecomm_reserv_richness)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.532 -1.853 -0.095  1.742  4.954 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 21.09296    1.91432  11.018 2.77e-08 ***
## shade        0.16474    0.09896   1.665    0.118    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.088 on 14 degrees of freedom
## Multiple R-squared:  0.1652, Adjusted R-squared:  0.1056 
## F-statistic: 2.771 on 1 and 14 DF,  p-value: 0.1182
Rich_reserv_hydro<-lm(interrich_n3~hydro,vegecomm_reserv_richness)
summary(Rich_reserv_hydro)
## 
## Call:
## lm(formula = interrich_n3 ~ hydro, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3095 -2.0680  0.4451  1.0590  5.1184 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  21.3598     1.1362  18.800 2.48e-11 ***
## hydro         0.6499     0.2252   2.885    0.012 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.676 on 14 degrees of freedom
## Multiple R-squared:  0.3729, Adjusted R-squared:  0.3281 
## F-statistic: 8.324 on 1 and 14 DF,  p-value: 0.01199
Rich_reserv_shrub<-lm(interrich_n3~shrub,vegecomm_reserv_richness)
summary(Rich_reserv_shrub)
## 
## Call:
## lm(formula = interrich_n3 ~ shrub, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5855 -2.2688  0.3838  1.5842  4.6146 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   22.147      1.055  20.996 5.56e-12 ***
## shrub          6.154      2.573   2.392   0.0314 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.848 on 14 degrees of freedom
## Multiple R-squared:  0.2901, Adjusted R-squared:  0.2394 
## F-statistic: 5.721 on 1 and 14 DF,  p-value: 0.03135
Rich_reserv_svp<-lm(interrich_n3~svp,vegecomm_reserv_richness)
summary(Rich_reserv_svp)
## 
## Call:
## lm(formula = interrich_n3 ~ svp, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8940 -2.3474  0.0939  2.4361  5.6107 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  23.2440     1.1032  21.069  5.3e-12 ***
## svp           0.9503     0.9240   1.028    0.321    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.259 on 14 degrees of freedom
## Multiple R-squared:  0.07025,    Adjusted R-squared:  0.003838 
## F-statistic: 1.058 on 1 and 14 DF,  p-value: 0.3212
Rich_reserv_ap<-lm(interrich_n3~ap,vegecomm_reserv_richness)
summary(Rich_reserv_ap)
## 
## Call:
## lm(formula = interrich_n3 ~ ap, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8181 -2.0705 -0.2523  1.5602  4.7442 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   21.675      1.532  14.152  1.1e-09 ***
## ap            25.248     14.352   1.759      0.1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058 on 14 degrees of freedom
## Multiple R-squared:  0.181,  Adjusted R-squared:  0.1225 
## F-statistic: 3.095 on 1 and 14 DF,  p-value: 0.1004
Rich_reserv_bryo<-lm(interrich_n3~bryo,vegecomm_reserv_richness)
summary(Rich_reserv_bryo)
## 
## Call:
## lm(formula = interrich_n3 ~ bryo, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7465 -1.9194  0.1346  1.4286  5.1210 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   23.096      1.047  22.065 2.83e-12 ***
## bryo          59.466     44.336   1.341    0.201    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.181 on 14 degrees of freedom
## Multiple R-squared:  0.1139, Adjusted R-squared:  0.05057 
## F-statistic: 1.799 on 1 and 14 DF,  p-value: 0.2012
Rich_reserv_sens_rel<-lm(interrich_n3~sens_rel,vegecomm_reserv_richness)
summary(Rich_reserv_sens_rel)
## 
## Call:
## lm(formula = interrich_n3 ~ sens_rel, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0560 -1.6362 -0.5409  1.4794  5.9762 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   22.042      1.382  15.953 2.25e-10 ***
## sens_rel      14.253      8.320   1.713    0.109    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.073 on 14 degrees of freedom
## Multiple R-squared:  0.1733, Adjusted R-squared:  0.1143 
## F-statistic: 2.935 on 1 and 14 DF,  p-value: 0.1087
Rich_reserv_tol_rel<-lm(interrich_n3~tol_rel,vegecomm_reserv_richness)
summary(Rich_reserv_tol_rel)
## 
## Call:
## lm(formula = interrich_n3 ~ tol_rel, data = vegecomm_reserv_richness)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.981 -1.519  0.258  1.307  3.891 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   28.411      1.599  17.765 5.33e-11 ***
## tol_rel      -10.939      3.623  -3.019  0.00919 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.63 on 14 degrees of freedom
## Multiple R-squared:  0.3943, Adjusted R-squared:  0.3511 
## F-statistic: 9.116 on 1 and 14 DF,  p-value: 0.009195
Rich_reserv_fqi<-lm(interrich_n3~fqi,vegecomm_reserv_richness)
summary(Rich_reserv_fqi)
## 
## Call:
## lm(formula = interrich_n3 ~ fqi, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1569 -1.5216  0.2207  1.6862  5.0045 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.5248     2.4319   7.617 2.41e-06 ***
## fqi           0.3606     0.1528   2.359   0.0334 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.859 on 14 degrees of freedom
## Multiple R-squared:  0.2845, Adjusted R-squared:  0.2334 
## F-statistic: 5.566 on 1 and 14 DF,  p-value: 0.03336
Rich_reserv_inv<-lm(interrich_n3~inv,vegecomm_reserv_richness)
summary(Rich_reserv_inv)
## 
## Call:
## lm(formula = interrich_n3 ~ inv, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4079 -2.2428  0.2233  1.6034  5.3265 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   23.758      0.951  24.983 5.17e-13 ***
## inv           27.430     49.518   0.554    0.588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.343 on 14 degrees of freedom
## Multiple R-squared:  0.02145,    Adjusted R-squared:  -0.04845 
## F-statistic: 0.3069 on 1 and 14 DF,  p-value: 0.5884
Rich_reserv_small_tree<-lm(interrich_n3~small_tree,vegecomm_reserv_richness)
summary(Rich_reserv_small_tree)
## 
## Call:
## lm(formula = interrich_n3 ~ small_tree, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.9406 -2.1955  0.3979  1.8062  4.2438 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    21.90       1.93  11.346 1.91e-08 ***
## small_tree     15.98      13.30   1.201     0.25    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.218 on 14 degrees of freedom
## Multiple R-squared:  0.09346,    Adjusted R-squared:  0.02871 
## F-statistic: 1.443 on 1 and 14 DF,  p-value: 0.2495
Rich_reserv_subcanopy<-lm(interrich_n3~subcanopy,vegecomm_reserv_richness)
summary(Rich_reserv_subcanopy)
## 
## Call:
## lm(formula = interrich_n3 ~ subcanopy, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1575 -1.2786  0.9301  1.9275  3.7641 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   22.639      1.269  17.838 5.04e-11 ***
## subcanopy     23.877     17.271   1.383    0.188    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.17 on 14 degrees of freedom
## Multiple R-squared:  0.1201, Adjusted R-squared:  0.05727 
## F-statistic: 1.911 on 1 and 14 DF,  p-value: 0.1885
Rich_reserv_canopy<-lm(interrich_n3~canopy,vegecomm_reserv_richness)
summary(Rich_reserv_canopy)
## 
## Call:
## lm(formula = interrich_n3 ~ canopy, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4648 -2.6753  0.5155  1.7726  5.3828 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   26.576      4.432   5.997 3.27e-05 ***
## canopy       -12.717     21.561  -0.590    0.565    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.338 on 14 degrees of freedom
## Multiple R-squared:  0.02424,    Adjusted R-squared:  -0.04545 
## F-statistic: 0.3479 on 1 and 14 DF,  p-value: 0.5647
Rich_reserv_null<-lm(interrich_n3~1,vegecomm_reserv_richness)
summary(Rich_reserv_null)
## 
## Call:
## lm(formula = interrich_n3 ~ 1, data = vegecomm_reserv_richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6592 -2.4610  0.2916  1.9212  5.1310 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  24.0092     0.8163   29.41 1.11e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.265 on 15 degrees of freedom
richmods_reserv<-list(Rich_reserv_ap,Rich_reserv_bryo,Rich_reserv_canopy,Rich_reserv_carex,Rich_reserv_cyper,Rich_reserv_dicot,Rich_reserv_fqi,Rich_reserv_hydro,Rich_reserv_inv,Rich_reserv_sens_rel,Rich_reserv_shade,Rich_reserv_shrub,Rich_reserv_small_tree,Rich_reserv_subcanopy,Rich_reserv_svp,Rich_reserv_tol_rel,Rich_reserv_vibifq,Rich_reserv_vibifqnotrees,Rich_reserv_null)

aictab(richmods_reserv)
## Warning in aictab.AIClm(richmods_reserv): 
## Model names have been supplied automatically in the table
## 
## Model selection based on AICc:
## 
##       K  AICc Delta_AICc AICcWt Cum.Wt     LL
## Mod18 3 80.88       0.00   0.36   0.36 -36.44
## Mod16 3 82.21       1.34   0.18   0.54 -37.11
## Mod8  3 82.77       1.89   0.14   0.68 -37.39
## Mod17 3 83.41       2.54   0.10   0.78 -37.71
## Mod12 3 84.76       3.88   0.05   0.83 -38.38
## Mod7  3 84.88       4.00   0.05   0.88 -38.44
## Mod1  3 87.04       6.16   0.02   0.89 -39.52
## Mod19 2 87.16       6.28   0.02   0.91 -41.12
## Mod10 3 87.19       6.31   0.02   0.92 -39.60
## Mod11 3 87.35       6.47   0.01   0.94 -39.67
## Mod6  3 87.97       7.09   0.01   0.95 -39.98
## Mod14 3 88.19       7.31   0.01   0.96 -40.10
## Mod2  3 88.30       7.42   0.01   0.97 -40.15
## Mod13 3 88.67       7.79   0.01   0.97 -40.33
## Mod4  3 88.73       7.85   0.01   0.98 -40.36
## Mod15 3 89.07       8.19   0.01   0.99 -40.54
## Mod5  3 89.17       8.29   0.01   0.99 -40.59
## Mod3  3 89.85       8.97   0.00   1.00 -40.92
## Mod9  3 89.89       9.01   0.00   1.00 -40.95

although our power is a bit limited above, it does appear that the Vegetative Index of Biotic Integrity - at a reservation-level - does appear to predict bird species richness (as long as we do not include woody plants in that list!)

we can also plot out richness for each site:

ggplot(aes(x=reservatio,y=reservspcount),data=vegecomm_reserv_richness)+geom_col()

ggplot(aes(x=reservatio,y=interrich_n3),data=vegecomm_reserv_richness)+geom_col()

In confirmation of our prior ordination plots, it now appears that Rocky River is one of the most species rich areas - followed by Hinckley, South Chagrin, Bedford, and Brecksville.

reservation shape features and species richness

While these vegetative associations are certainly important, it is often helpful to also think about how the shape of our reservations influences the species found within them. There is a rich body of literature on the effects of landscape context and configuration on

vegecomm_reserv_richness
##    reservatio    vibifq vibifqnotrees     carex     cyper     dicot     shade     shrub     hydro       svp         ap       fqi        bryo   hydro_rel hydro_rel_notrees    sens_rel sens_rel_notree   tol_rel tol_rel_notree          inv small_tree  subcanopy    canopy numsurveys reservspcount
## BC         BC 34.601035     25.313664 3.8158996 4.1087866 16.669456 15.774059 0.3807531 5.6276151 0.0000000 0.11874312 14.284368 0.004958120 0.023797808       0.028223440 0.089812897     0.045449261 0.4891651      0.5660586 0.0295276552 0.17029270 0.01723615 0.1166671         18            48
## BE         BE 39.504613     29.726028 3.0626566 3.1102757 21.440267 20.778613 0.2063492 2.8671679 1.4152047 0.14248544 15.672360 0.014829724 0.047829388       0.078499350 0.164717698     0.070814554 0.3887067      0.4714974 0.0051759616 0.11508660 0.07952218 0.1801468        100            73
## BK         BK 40.504890     33.317087 0.6134454 0.6134454  9.773109  6.680672 0.0000000 0.7731092 0.0000000 0.18309632 13.396361 0.081362435 0.000000000       0.000000000 0.366516404     0.310637023 0.1836977      0.1760264 0.0000000000 0.09542484 0.11412427 0.1443458          9            45
## BR         BR 55.285843     41.682545 5.2304060 5.4381492 28.100094 30.229462 0.4277620 4.4711992 2.7648725 0.10423918 21.496159 0.022946191 0.035888066       0.070085339 0.166597534     0.103523925 0.3088704      0.3869337 0.0564445807 0.18449116 0.05321523 0.2079372         82            70
## BW         BW 47.464208     35.190650 3.4246285 3.5095541 18.643312 18.087049 0.4394904 2.7940552 1.9341826 0.07325213 19.131091 0.016114563 0.049254626       0.106742227 0.160002194     0.044342921 0.2891663      0.3490757 0.0000000000 0.15571638 0.08904453 0.1886476         41            41
## EC         EC 31.350954     11.944005 0.0000000 0.0000000 19.000000 10.000000 0.0000000 3.0000000 0.0000000 0.00000000 12.247449 0.012371644 0.000000000       0.000000000 0.024743288     0.000000000 0.4353994      0.6448408 0.0000000000 0.16216216 0.14737364 0.2303164          3            20
## GP         GP 27.946924     10.651350 2.0000000 2.0000000 23.000000 16.000000 0.0000000 7.0000000 0.0000000 0.00000000 12.809880 0.009602612 0.016260423       0.025736505 0.177936399     0.000000000 0.4184018      0.6011348 0.0096026120 0.06153846 0.00000000 0.2061372          3            24
## HI         HI 56.777233     46.984497 4.2806748 4.3550613 38.047546 34.980061 0.7730061 9.4685583 1.6541411 0.08149579 23.192002 0.009410177 0.088963035       0.127416913 0.107577765     0.055443524 0.2894462      0.3644430 0.0003946167 0.14618568 0.09429889 0.1791128         92            78
## HU         HU 24.435647      6.752597 3.0000000 4.0000000 30.000000 24.000000 0.0000000 0.0000000 0.0000000 0.10526316 14.293825 0.004343231 0.000000000       0.000000000 0.078178159     0.000000000 0.6169125      0.8298138 0.0000000000 0.04727273 0.00000000 0.2077735          7            27
## MS         MS 46.355466     29.687439 3.5869738 3.7172359 20.918983 19.845909 0.3494837 4.3288324 0.7648928 0.08450301 17.252813 0.011838657 0.048039227       0.123442813 0.165213079     0.032381386 0.2526545      0.3385182 0.0414410260 0.14354510 0.05094113 0.1931117        104            63
## NC         NC 53.965315     38.560833 2.3760684 2.4658120 23.248718 20.439316 0.8153846 7.8316239 1.3606838 0.10368316 19.929270 0.014654907 0.112169719       0.209235627 0.227754177     0.127802382 0.2922618      0.3813379 0.0002100554 0.12328842 0.06749561 0.2470918         96            66
## OE         OE 13.813003     14.088671 3.6745283 3.6745283 17.641509 14.811321 0.0000000 0.1839623 0.8207547 0.16659520 10.636789 0.002949977 0.001007902       0.001672384 0.005627719     0.006567877 0.5388112      0.6840616 0.0016111062 0.05005438 0.00000000 0.1768963         18            43
## RR         RR 39.382635     27.503422 1.4131868 1.4717949 23.423443 18.774359 0.3897436 4.0959707 0.3003663 0.11351880 15.560389 0.015515483 0.084112185       0.123459302 0.078730158     0.025448674 0.2574988      0.3879014 0.0020324549 0.19507519 0.11464307 0.2216719         93            80
## SC         SC 39.740571     28.655854 3.1097209 3.3666987 22.470645 18.662175 0.6496631 8.8960539 1.8671800 0.07591598 16.413909 0.010576548 0.072576642       0.128944299 0.147712040     0.073316453 0.4357619      0.5385634 0.0001414818 0.08700562 0.06758911 0.2247417         76            76
## WA         WA  1.137611      1.026171 0.0000000 0.0000000  4.000000  3.000000 0.0000000 0.0000000 0.0000000 0.00000000  2.474874 0.000000000 0.000000000       0.000000000 0.000000000     0.000000000 0.9516582      0.9604965 0.0000000000 0.08695652 0.00000000 0.2171573          6            23
## WC         WC 39.518844     23.929770 3.2511628 3.6604651 13.069767 11.162791 0.4093023 3.8883721 0.0000000 0.12656059 14.560841 0.014097267 0.126317347       0.209907312 0.247360772     0.065370632 0.2894291      0.2196096 0.0000000000 0.28662562 0.02284111 0.2881646         17            48
##    interrich_n3
## BC     25.65319
## BE     25.03640
## BK     26.76190
## BR     25.61177
## BW     20.26632
## EC     20.00000
## GP     24.00000
## HI     28.47655
## HU     19.51429
## MS     23.31928
## NC     23.90454
## OE     21.97549
## RR     29.14014
## SC     27.53519
## WA     18.35000
## WC     24.60147
structuralconnectivity<-read.csv("Reservations_LandMetrics_bbs.csv",header=T)
vegecomm_reserv_richness$reservatio
##  [1] "BC" "BE" "BK" "BR" "BW" "EC" "GP" "HI" "HU" "MS" "NC" "OE" "RR" "SC" "WA" "WC"
structuralconnectivity$resID
##  [1] "BC" "BE" "BK" "BR" "BW" "EC" "GP" "HI" "HU" "MS" "NC" "OE" "RR" "SC" "WA" "WC"
vegecomm_reserv_richness_struct<-cbind(vegecomm_reserv_richness,structuralconnectivity)

Rich_reserv_acres<-lm(interrich_n3~acres,vegecomm_reserv_richness_struct)
summary(Rich_reserv_acres)
## 
## Call:
## lm(formula = interrich_n3 ~ acres, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7093 -2.6366 -0.5477  2.3302  3.5818 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 21.977011   1.013861  21.677  3.6e-12 ***
## acres        0.001382   0.000509   2.715   0.0168 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.735 on 14 degrees of freedom
## Multiple R-squared:  0.3449, Adjusted R-squared:  0.2981 
## F-statistic: 7.371 on 1 and 14 DF,  p-value: 0.01676
Rich_reserv_ai<-lm(interrich_n3~ai,vegecomm_reserv_richness_struct)
summary(Rich_reserv_ai)
## 
## Call:
## lm(formula = interrich_n3 ~ ai, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1392 -2.0183  0.6307  2.3067  4.5389 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   4.4278    14.9619   0.296    0.772
## ai            0.2097     0.1600   1.311    0.211
## 
## Residual standard error: 3.19 on 14 degrees of freedom
## Multiple R-squared:  0.1093, Adjusted R-squared:  0.04566 
## F-statistic: 1.718 on 1 and 14 DF,  p-value: 0.2111
Rich_reserv_areamn<-lm(interrich_n3~area_mn,vegecomm_reserv_richness_struct)
summary(Rich_reserv_areamn)
## 
## Call:
## lm(formula = interrich_n3 ~ area_mn, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2080 -1.9539  0.2301  2.4859  4.5089 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.251e+01  9.114e-01  24.697 6.06e-13 ***
## area_mn     1.297e-03  5.071e-04   2.558   0.0228 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.79 on 14 degrees of freedom
## Multiple R-squared:  0.3184, Adjusted R-squared:  0.2698 
## F-statistic: 6.541 on 1 and 14 DF,  p-value: 0.02278
Rich_reserv_ca<-lm(interrich_n3~ca,vegecomm_reserv_richness_struct)
summary(Rich_reserv_ca)
## 
## Call:
## lm(formula = interrich_n3 ~ ca, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7135 -2.6368 -0.5517  2.3337  3.5862 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.198e+01  1.014e+00  21.686 3.58e-12 ***
## ca          3.169e-04  1.168e-04   2.714   0.0168 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.736 on 14 degrees of freedom
## Multiple R-squared:  0.3447, Adjusted R-squared:  0.2979 
## F-statistic: 7.363 on 1 and 14 DF,  p-value: 0.01681
Rich_reserv_cai_mn<-lm(interrich_n3~cai_mn,vegecomm_reserv_richness_struct)
summary(Rich_reserv_cai_mn)
## 
## Call:
## lm(formula = interrich_n3 ~ cai_mn, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8640 -1.2744  0.0397  1.4097  4.9520 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 22.01786    1.69266  13.008  3.3e-09 ***
## cai_mn       0.05042    0.03783   1.333    0.204    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.184 on 14 degrees of freedom
## Multiple R-squared:  0.1126, Adjusted R-squared:  0.04924 
## F-statistic: 1.777 on 1 and 14 DF,  p-value: 0.2038
Rich_reserv_circle_mn<-lm(interrich_n3~circle_mn,vegecomm_reserv_richness_struct)
summary(Rich_reserv_circle_mn)
## 
## Call:
## lm(formula = interrich_n3 ~ circle_mn, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.2121 -1.0889 -0.1669  2.3594  4.4504 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   13.241      5.059   2.617   0.0203 *
## circle_mn     16.498      7.670   2.151   0.0494 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.93 on 14 degrees of freedom
## Multiple R-squared:  0.2484, Adjusted R-squared:  0.1947 
## F-statistic: 4.627 on 1 and 14 DF,  p-value: 0.04942
Rich_reserv_clumpy<-lm(interrich_n3~clumpy,vegecomm_reserv_richness_struct)
summary(Rich_reserv_clumpy)
## 
## Call:
## lm(formula = interrich_n3 ~ clumpy, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5004 -1.3690  0.5353  1.5123  3.6057 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   131.88      32.28   4.085  0.00111 **
## clumpy       -105.17      31.47  -3.342  0.00484 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.521 on 14 degrees of freedom
## Multiple R-squared:  0.4438, Adjusted R-squared:  0.404 
## F-statistic: 11.17 on 1 and 14 DF,  p-value: 0.004841
Rich_reserv_cohesion<-lm(interrich_n3~cohesion,vegecomm_reserv_richness_struct)
summary(Rich_reserv_cohesion)
## 
## Call:
## lm(formula = interrich_n3 ~ cohesion, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7345 -1.7224  0.3427  1.9993  3.2565 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -64.3888    30.6171  -2.103   0.0540 .
## cohesion      0.9087     0.3146   2.888   0.0119 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.675 on 14 degrees of freedom
## Multiple R-squared:  0.3733, Adjusted R-squared:  0.3286 
## F-statistic:  8.34 on 1 and 14 DF,  p-value: 0.01192
Rich_reserv_contig_mn<-lm(interrich_n3~contig_mn,vegecomm_reserv_richness_struct)
summary(Rich_reserv_contig_mn)
## 
## Call:
## lm(formula = interrich_n3 ~ contig_mn, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8000 -1.5313 -0.0424  1.6983  5.5598 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   20.330      2.339   8.693 5.14e-07 ***
## contig_mn      6.071      3.643   1.667    0.118    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.087 on 14 degrees of freedom
## Multiple R-squared:  0.1655, Adjusted R-squared:  0.1059 
## F-statistic: 2.777 on 1 and 14 DF,  p-value: 0.1178
Rich_reserv_core_mn<-lm(interrich_n3~core_mn,vegecomm_reserv_richness_struct)
summary(Rich_reserv_core_mn)
## 
## Call:
## lm(formula = interrich_n3 ~ core_mn, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2463 -1.9870  0.2502  2.3776  4.5108 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.257e+01  8.966e-01   25.17 4.67e-13 ***
## core_mn     1.470e-03  5.742e-04    2.56   0.0227 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.789 on 14 degrees of freedom
## Multiple R-squared:  0.3188, Adjusted R-squared:  0.2702 
## F-statistic: 6.553 on 1 and 14 DF,  p-value: 0.02268
Rich_reserv_cpland<-lm(interrich_n3~cpland,vegecomm_reserv_richness_struct)
summary(Rich_reserv_cpland)
## 
## Call:
## lm(formula = interrich_n3 ~ cpland, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8173 -2.7253 -0.4908  2.2622  3.6100 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.1180     0.9964  22.198 2.61e-12 ***
## cpland        0.3716     0.1411   2.634   0.0196 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.764 on 14 degrees of freedom
## Multiple R-squared:  0.3313, Adjusted R-squared:  0.2836 
## F-statistic: 6.938 on 1 and 14 DF,  p-value: 0.01963
Rich_reserv_dcad<-lm(interrich_n3~dcad,vegecomm_reserv_richness_struct)
summary(Rich_reserv_dcad)
## 
## Call:
## lm(formula = interrich_n3 ~ dcad, data = vegecomm_reserv_richness_struct)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4125 -2.8268  0.1998  1.4137  5.5563 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   22.645      1.201   18.86 2.39e-11 ***
## dcad         124.637     83.107    1.50    0.156    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.137 on 14 degrees of freedom
## Multiple R-squared:  0.1384, Adjusted R-squared:  0.07688 
## F-statistic: 2.249 on 1 and 14 DF,  p-value: 0.1559

multi-species occupancy modelling

MSOM - running different models

basic null model

bbs_occobject_null.occ.formula <- ~1
bbs_occobject_null.det.formula <- ~ 1
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.null <- msPGOcc(occ.formula = bbs_occobject_null.occ.formula, 
                  det.formula = bbs_occobject_null.det.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.null)
ppc.ms.null.out <- ppcOcc(out.ms.null, 'chi-squared', group = 1)
summary(ppc.ms.null.out)

waicOcc(out.ms.null)

building up complexity of detection model

bbs_occobject_null.occ.formula <- ~1
bbs_occobject_day.det.formula <- ~ scale(day)
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.detect.day <- msPGOcc(occ.formula = bbs_occobject_null.occ.formula, 
                  det.formula = bbs_occobject_day.det.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.detect.day)
ppc.ms.detect.day.out <- ppcOcc(out.ms.detect.day, 'chi-squared', group = 1)
summary(out.ms.detect.day)

waicOcc(out.ms.detect.day)

day and time-of-day

bbs_occobject_null.occ.formula <- ~1
bbs_occobject_daytod.det.formula <- ~ scale(day)+scale(tod)
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.detect.daytod <- msPGOcc(occ.formula = bbs_occobject_null.occ.formula, 
                  det.formula = bbs_occobject_daytod.det.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.detect.daytod)
ppc.ms.detect.day.out <- ppcOcc(out.ms.detect.daytod, 'chi-squared', group = 1)
summary(out.ms.detect.daytod)

waicOcc(out.ms.detect.daytod)

day, tod, quadratic day

bbs_occobject_null.occ.formula <- ~1
bbs_occobject_daytodquad.det.formula <- ~ scale(day)+scale(tod) +I(scale(day)^2)
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.detect.daytodquad<- msPGOcc(occ.formula = bbs_occobject_null.occ.formula, 
                  det.formula = bbs_occobject_daytodquad.det.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.detect.daytodquad)
ppc.ms.detect.day.out <- ppcOcc(out.ms.detect.daytodquad, 'chi-squared', group = 1)
summary(out.ms.detect.daytodquad)

waicOcc(out.ms.detect.daytodquad)

The best-supported model is within 2 WAIC of a simpler detection model - one that just includes day. So we will proceed with just using a linear term for day in our detectability submodels.

Now we will run several different occupancy submodels, once again evaluating fit for progressively more complicated models.

occupancy model with vibi-fq

bbs_occobject_vibifq.occ.formula <- ~ scale(vibifq)
bbs_occobject.det.best.formula <- ~ scale(day) 
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.vibifq <- msPGOcc(occ.formula = bbs_occobject_vibifq.occ.formula, 
                  det.formula = bbs_occobject.det.best.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.vibifq)
ppc.ms.vibifq.out <- ppcOcc(out.ms.vibifq, 'chi-squared', group = 1)
summary(ppc.ms.vibifq.out)


out.ms.bestdec.small <- msPGOcc(occ.formula = ~1, 
                       det.formula = bbs_occobject.det.best.formula, 
                       data = bbs_occobject, 
                       inits = ms.inits, 
                       n.samples = 30000, 
                       priors = ms.priors, 
                       n.omp.threads = 1, 
                       verbose = TRUE, 
                       n.report = 6000, 
                       n.burn = 10000,
                       n.thin = 50, 
                       n.chains = 3)


waicOcc(out.ms.bestdec.small)
waicOcc(out.ms.vibifq)

occupancy model with vibi-fq-notrees

bbs_occobject_vibifqnotrees.occ.formula <- ~ scale(vibifqnotrees)
bbs_occobject.det.best.formula <- ~ scale(day) 
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.vibifqnotrees <- msPGOcc(occ.formula = bbs_occobject_vibifqnotrees.occ.formula, 
                  det.formula = bbs_occobject.det.best.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.vibifqnotrees)
ppc.ms.vibifqnotrees.out <- ppcOcc(out.ms.vibifqnotrees, 'chi-squared', group = 1)
summary(ppc.ms.vibifqnotrees.out)

waicOcc(out.ms.vibifqnotrees)

occupancy model with all rda-supported covariates

bbs_occobject_fullrda.occ.formula <- ~ scale(tol_rel) + scale(svp) + scale(dicot) + scale(cyper) + scale(small_tree)
bbs_occobject.det.best.formula <- ~ scale(day) 
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.fullrda <- msPGOcc(occ.formula = bbs_occobject_fullrda.occ.formula, 
                  det.formula = bbs_occobject.det.best.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.fullrda)
ppc.ms.fullrda.out <- ppcOcc(out.ms.fullrda, 'chi-squared', group = 1)
summary(ppc.ms.fullrda.out)

waicOcc(out.ms.fullrda)

occupancy model with just tolerance/sensitivity

bbs_occobject_tolerance.occ.formula <- ~ scale(tol_rel)
bbs_occobject.det.best.formula <- ~ scale(day) 
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.tolerance<- msPGOcc(occ.formula = bbs_occobject_tolerance.occ.formula, 
                  det.formula = bbs_occobject.det.best.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.tolerance)
ppc.ms.tolerance.out <- ppcOcc(out.ms.tolerance, 'chi-squared', group = 1)
summary(ppc.ms.tolerance.out)

waicOcc(out.ms.tolerance)

global model

bbs_occobject_global.occ.formula <- ~ scale(tol_rel) + scale(svp) + scale(dicot) + scale(cyper) + scale(small_tree)+scale(shrub)+scale(hydro)+scale(subcanopy)+scale(canopy)
bbs_occobject.det.best.formula <- ~ scale(day) 
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.global <- msPGOcc(occ.formula = bbs_occobject_global.occ.formula, 
                  det.formula = bbs_occobject.det.best.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.global)
ppc.ms.global.out <- ppcOcc(out.ms.global, 'chi-squared', group = 1)
summary(ppc.ms.global.out)

waicOcc(out.ms.global)

global model - all non-correlated covariates (after removing some redundancy)

bbs_occobject_global.occ.formula <- ~ scale(tol_rel) + scale(svp) + scale(dicot) + scale(cyper) + scale(small_tree)+scale(shrub)+scale(hydro)+scale(subcanopy)+scale(canopy)
bbs_occobject.det.best.formula <- ~ scale(day) 
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.global <- msPGOcc(occ.formula = bbs_occobject_global.occ.formula, 
                  det.formula = bbs_occobject.det.best.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.global)
ppc.ms.global.out <- ppcOcc(out.ms.global, 'chi-squared', group = 1)
summary(ppc.ms.global.out)

waicOcc(out.ms.global)

just FQAI

bbs_occobject_fqi.occ.formula <- ~ scale(fqi)
bbs_occobject.det.best.formula <- ~ scale(day) 
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0, 
                 beta.comm = 0, 
                 beta = 0, 
                 alpha = 0,
                 tau.sq.beta = 1, 
                 tau.sq.alpha = 1, 
                 z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
                  alpha.comm.normal = list(mean = 0, var = 2.72), 
                  tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                  tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
out.ms.fqi <- msPGOcc(occ.formula = bbs_occobject_fqi.occ.formula, 
                  det.formula = bbs_occobject.det.best.formula, 
                  data = bbs_occobject, 
                  inits = ms.inits, 
                  n.samples = 30000, 
                  priors = ms.priors, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 6000, 
                  n.burn = 10000,
                  n.thin = 50, 
                  n.chains = 3)


summary(out.ms.fqi)
ppc.ms.fqi.out <- ppcOcc(out.ms.fqi, 'chi-squared', group = 1)
summary(ppc.ms.fqi.out)

waicOcc(out.ms.fqi)

model selection

we can also print out a table of Widely Applicable Information Criterion (wAIC) values. The lower this is, the better the model fit.

AICc_table<-tibble(.rows = 10)

AICc_table$model<-c("Null","p ~ day", "p ~ day + time-of-day", "p ~ day^2 + day +tod", "p ~ day, psi ~ vibi-fq", "p ~ day, psi ~ vibi-fq-notree","p ~ day, psi ~ tolerant + seedless vascular + dicots + cyperaceae + small trees","p~day,psi~tolerant","p~day,psi~tolerant + seedless vascular + dicot + cyperaceae + small trees +shrubs+hydro+subcanopy+canopy","p~day,psi~fqi")

AICc_table$wAIC<-c(waicOcc(out.ms.null)[[3]],waicOcc(out.ms.detect.day)[[3]],waicOcc(out.ms.detect.daytod)[[3]],
                   waicOcc(out.ms.detect.daytodquad)[[3]],waicOcc(out.ms.vibifq)[[3]],waicOcc(out.ms.vibifqnotrees)[[3]],
                   waicOcc(out.ms.fullrda)[[3]],waicOcc(out.ms.tolerance)[[3]],waicOcc(out.ms.global)[[3]],waicOcc(out.ms.fqi))
##   X                                                                                                    model     wAIC
## 9 9 p~day,psi~tolerant + seedless vascular + dicot + cyperaceae + small trees +shrubs+hydro+subcanopy+canopy 23165.73
## 7 7                          p ~ day, psi ~ tolerant + seedless vascular + dicots + cyperaceae + small trees 23185.00
## 5 5                                                                                   p ~ day, psi ~ vibi-fq 23190.57
## 8 8                                                                                       p~day,psi~tolerant 23323.25
## 6 6                                                                            p ~ day, psi ~ vibi-fq-notree 23372.44
## 4 4                                                                                     p ~ day^2 + day +tod 23755.71
## 2 2                                                                                                  p ~ day 23756.93
## 3 3                                                                                    p ~ day + time-of-day 23769.78
## 1 1                                                                                                     Null 23840.32

So, the top model is the global, followed by the fullrda model.

Now, we plot based on this top model.

plotting using an RShiny app

covariate plotting

#options(max.print=1000000)
#sink("outmsglobal_90percent.txt")
#print(summary(out.ms.global,level="species",quantiles=c(0.05,0.5,0.95)))
#sink() 

covars<-read.csv("outmsglobal_covar.csv",header=T)

covars_tolerance<-subset(covars,covars$Covariate=="tol_rel")

covars_tolerance<-covars_tolerance[order(covars_tolerance$Mean),]

spp_tolorder<-factor(covars_tolerance$Species,levels=covars_tolerance$Species)

covars_tolerance$spporder<-spp_tolorder

ggplot(aes(x=Mean,y=spp_tolorder),data=covars_tolerance)+geom_point()+geom_errorbarh(aes(xmin=X5.0CI,xmax=X95.0CI))+geom_vline(xintercept=0.1652)+theme_classic(base_size=16)

covars_svp<-subset(covars,covars$Covariate=="svp")

covars_svp<-covars_svp[order(covars_svp$Mean),]


spp_svporder<-factor(covars_svp$Species,levels=covars_svp$Species)

covars_svp$spporder<-spp_svporder

ggplot(aes(x=Mean,y=spp_svporder),data=covars_svp)+geom_point()+geom_errorbarh(aes(xmin=X5.0CI,xmax=X95.0CI))+geom_vline(xintercept=-0.0594)+theme_classic(base_size=16)

covars_dicot<-subset(covars,covars$Covariate=="dicot")

covars_dicot<-covars_dicot[order(covars_dicot$Mean),]


spp_dicotorder<-factor(covars_dicot$Species,levels=covars_dicot$Species)

covars_dicot$spporder<-spp_dicotorder

ggplot(aes(x=Mean,y=spp_dicotorder),data=covars_dicot)+geom_point()+geom_errorbarh(aes(xmin=X5.0CI,xmax=X95.0CI))+geom_vline(xintercept=0.2461)+theme_classic(base_size=16)

covars_cyper<-subset(covars,covars$Covariate=="cyper")

covars_cyper<-covars_cyper[order(covars_cyper$Mean),]


spp_cyperorder<-factor(covars_cyper$Species,levels=covars_cyper$Species)

covars_cyper$spporder<-spp_cyperorder

ggplot(aes(x=Mean,y=spp_cyperorder),data=covars_cyper)+geom_point()+geom_errorbarh(aes(xmin=X5.0CI,xmax=X95.0CI))+geom_vline(xintercept=-0.1375)+theme_classic(base_size=16)

covars_smalltree<-subset(covars,covars$Covariate=="small_tree")

covars_smalltree<-covars_smalltree[order(covars_smalltree$Mean),]


spp_treeorder<-factor(covars_smalltree$Species,levels=covars_smalltree$Species)

covars_smalltree$spporder<-spp_treeorder

ggplot(aes(x=Mean,y=spp_treeorder),data=covars_smalltree)+geom_point()+geom_errorbarh(aes(xmin=X5.0CI,xmax=X95.0CI))+geom_vline(xintercept=0.0672)+theme_classic(base_size=16)

covars_shrub<-subset(covars,covars$Covariate=="shrub")

covars_shrub<-covars_shrub[order(covars_shrub$Mean),]


spp_shruborder<-factor(covars_shrub$Species,levels=covars_shrub$Species)

covars_shrub$spporder<-spp_shruborder

ggplot(aes(x=Mean,y=spp_shruborder),data=covars_shrub)+geom_point()+geom_errorbarh(aes(xmin=X5.0CI,xmax=X95.0CI))+geom_vline(xintercept=-0.2181)+theme_classic(base_size=16)

covars_hydro<-subset(covars,covars$Covariate=="hydro")

covars_hydro<-covars_hydro[order(covars_hydro$Mean),]

spp_hydroorder<-factor(covars_hydro$Species,levels=covars_hydro$Species)

covars_hydro$spporder<-spp_hydroorder

ggplot(aes(x=Mean,y=spp_hydroorder),data=covars_hydro)+geom_point()+geom_errorbarh(aes(xmin=X5.0CI,xmax=X95.0CI))+geom_vline(xintercept=0.2068)+theme_classic(base_size=16)

covars_subcanopy<-subset(covars,covars$Covariate=="subcanopy")

covars_subcanopy<-covars_subcanopy[order(covars_subcanopy$Mean),]

spp_subcanopyorder<-factor(covars_subcanopy$Species,levels=covars_subcanopy$Species)

covars_subcanopy$spporder<-spp_subcanopyorder

ggplot(aes(x=Mean,y=spp_subcanopyorder),data=covars_subcanopy)+geom_point()+geom_errorbarh(aes(xmin=X5.0CI,xmax=X95.0CI))+geom_vline(xintercept=-0.0393)+theme_classic(base_size=16)

covars_canopy<-subset(covars,covars$Covariate=="canopy")

covars_canopy<-covars_canopy[order(covars_canopy$Mean),]

spp_canopyorder<-factor(covars_canopy$Species,levels=covars_canopy$Species)

covars_canopy$spporder<-spp_canopyorder

ggplot(aes(x=Mean,y=spp_canopyorder),data=covars_canopy)+geom_point()+geom_errorbarh(aes(xmin=X5.0CI,xmax=X95.0CI))+geom_vline(xintercept=-0.0702)+theme_classic(base_size=16)